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EVALUATION  OF  THE  GENERALIZED  QM  FUNCTION 
FOR  COMPLEX  ARGUMENTS 


INTRODUCTION 


The  generalized  QM  function  is  defined  as  (reference  1  and  reference  2,  section  VII, 
equation  2.18) 


WJ 

QM(a,b)  =  J  dxx 


\aj 


exp 


f  x  2+a2^ 


Im-M  *)• 


(1) 


Typically,  Mis  an  integer,  while  parameters  a  and  b  are  real  and  nonnegative.  Numerous 
integrals  involving  the  Qu(a,b )  function  are  available  in  references  3  and  4.  Also,  a 
generalization  to  a  general  power  of  ( x/a ),  as  well  as  a  general  order  of  the  modified  Bessel 
function  I v,  was  made  in  references  3  and  5.  Numerical  evaluation  of  QM  ( a,b )  for  large  values 
of  arguments  a  and/or  b  is  often  plagued  by  underflow  or  overflow,  due  to  the  presence  of 
exponentials  with  large  positive  or  negative  real  arguments.  Furthermore,  in  recent  applications 
to  active  sonar  processing  (references  6  and  7),  the  need  for  evaluation  of  QM  ( a,b )  with 
complex  arguments  a  and  b  was  encountered;  in  particular,  the  joint  probability  density  function 
of  the  M  - 1  largest  envelope-squared  values  from  a  matched  filter  and  the  sum  of  the  remainder 
required  exactly  this  quantity  QM(a,b).  However,  the  usual  series  expansions  for  QM{a,b )  can 
now  contain  powers  of  large  complex  numbers,  thereby  destroying  any  significance  in  the  final 
result,  due  to  cancellations  of  large  positive  and  negative  values  of  the  real  and  imaginary 
components  of  the  terms  in  the  series. 

To  construct  a  viable  routine  for  evaluation  of  QM  ( a,b )  with  arbitrary  complex  arguments, 
a  study  of  several  possibilities  was  undertaken.  The  study  begins  with  an  investigation  of  the 
basic  statistics  of  non-central  chi-squared  random  variables  (RVs),  and  progresses  into  the  realm 
of  probability  density  functions  (PDFs),  cumulative  distribution  functions  (CDFs),  and 
exceedance  distribution  functions  (EDFs),  where  the  QM(a,b)  function  is  encountered. 
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STATISTICS  OF  NON-CENTRAL  CHI-SQUARED  RANDOM  VARIABLES 


Let  {g(&)}  be  a  set  of  independent  Gaussian  RVs  with  means  { ju(k )}  and  standard 
deviations  {<j (k)},  for  k  =  1 :  K.  Then,  the  RVs 

y (k)  =  g (k)2  for  k  =  \:K 

have  the  moment- generating  functions  (MGFs) 


=  £{exp(Ay(A:))}  =  £{exp(;ig(£)2)}  =  f du  exp (Au2)  -=J-  —  exp 

J  a(k)  ^  2cr(k) 


=  (\  —  2a(k)2  X)~V2  exp 


1-2  a(ky  A 


2*rm  exp  1  for  Re(/t)<- — 7—7,  k  =  l  :K. 


2  cr(ky 


The  non-central  chi-squared  RV  of  interest  here  is  given  by  the  sum 


s = E  y(*) = E 


The  MGF  of  sum  s  follows  from  equation  (3)  and  the  independence  of  RVs  {g (&)}  as 


M,W  =  t\0-WW'n  «pft,  for  Re(A) <  — - f— -5-. 

l  7^1  1-2 a(k)  A)  2  max { cr{k )  } 


SPECIAL  CASE 

Let  the  standard  deviations  of  all  the  original  Gaussian  RVs  {g(£)}  be  equal;  that  is, 
a(k)  =  cr  for  k  =  1 :  K. 


Then,  the  MGF  in  equation  (5)  simplifies  to 


AW  =  (1  -2cr2AyK/2  exp  for  Re(A)  <  — ^ 

1  2.  (7  /t  2>Cf 


where 


m2  =  J^ju(ky 


2 


Observe  that  the  statistics  of  non-central  chi-squared  RV  s  depend  only  on  the  sum  m2  of  the 
squares  of  the  means,  {//(&)},  of  the  individual  Gaussian  RVs  {g(k)\  when  equation  (6)  is  true. 

The  PDF  of  RV  s ,  corresponding  to  MGF  (7),  is 


1  (J7A  f  u  +  m 2N1  (  m  r-\ 

ps(u)  =  p(u;K,m,o)  =  - — -  - —  exp  — -  IM_X  ~4u  for  u  >  0,  m>  0,(9) 

2cj  m  2(7  \o 


where 


M  =  m=  Yj  MkY 

2*  v  i=\ 


M  need  not  be  an  integer.  A  better  numerical  alternative  to  equation  (9)  is 


1  f  f  (■ yfu-m )2^  f  m^fu)  fmy/u) 

pM)  =  ~ T  -  exP - -  exp - T~  !M-i  — —  for  u>0’  m> °>  (n) 

2(7  ml  2<J  (7  a 


because  the  well-behaved  combination 


exp(-x)  =  besseli(M  - l,x,l) 


is  directly  available  as  a  MATLAB  function.  For  m  =  0,  equation  (11)  reduces  to 


Ps(U)  = 


r  (M)  (2cr  2y 


exp - -  for  u  >  0. 

2a2 


The  EDF  for  RV  s  in  equation  (4)  is,  upon  use  of  equations  (9)  and  (1), 


es(u)  =  Prob(s  >  u)  =  f  dt  ps(t)  =  QM  — for  u  >  0. 

J  (7  0 


The  y'-th  cumulant  of  RV  s  is 


(  M  m2  ^ 

Z.U)  =  jl  — +  (2a2y  for  7  =  1,2,..., 

V  J  2cr 


while  the  v  -th  moment  of  RV  s  is 
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.2,vr(M  +  v)  _  m 


E{sv}  =  (2cr2y  y- 

T{M)  1  l  2<r2 


in  terms  of  a  confluent  hypergeometric  function.  For  v  equal  to  an  integer  /,  this  result 
simplifies  to 


£{•'}  =  J!  (2ct2)' 

2  O’ 


where  the  latter  function  is  a  Laguerre  polynomial. 

An  alternative  RV  of  interest  is  the  square  root  of  the  sum  of  squares  of  the  Gaussian  RVs, 
namely,  non-central  chi-variate 

e  =  ^  =  ,Zg»2-  <18) 

V  k= 1 

For  the  case  of  equal  Gaussian  standard  deviations,  equation  (6),  the  PDF  of  RV  e  follows  from 
equation  (9)  as 


f  \M-l  Cl  l\  /  \ 

,  .  „  ,  u  u  \  u  +m  .  mu]  r  _  „ 

/?,(«)  =  2 «/?,(«  )  =  —  —  exp - — j—  V,  — r  for  u  >  0,  m  >  0. 

a  )  l  2cr  V  cr" 


The  corresponding  EDF  and  CDF  are,  respectively, 


.  .  _  m  u  ,  .  ,  /  \  i  si  t  m  ti  \  m  u  .... 

et(u)  =  QM  — ,  c,(i/)  =  l-ee(H)  =  l-0j*  — =  qM  — forn>0,  (20) 

VO’O’y  VO’O’y  VO’O’y 


where  qM  is  the  complementary  function.  A  better  numerical  alternative  to  equation  (19)  is 


«  f h)  f  (u-m)2  |  (  mu]  ( mu] 

pe{u)  =  —  —  exp - -  —  exp  -—  IMA  —  for  «  >  0,  w>0; 

cr  Im  2cr  V  cr  1  \a  J 


see  equation  (12).  For  m  =  0,  equation  (21)  reduces  to 


2u2M~]  u2  \ 

Pe(u)  = - -z—rr  exp - —  for  u  >  0. 

T(M)  (2cr2)M  \  2a2) 
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DETERMINATION  OF  THE  MIXTURE  FUNCTION 


A  mixture  of  a  CDF  and  an  MGF  was  defined  in  reference  6,  equation  (3),  namely, 

U 

c(u,A)  =  J  dx  p(x)  exp(A  x),  (23) 

-oo 

where  p(x )  is  a  PDF.  That  is,  c(u, 0)  is  the  usual  CDF,  while  c(+co,A)  is  the  usual  MGF. 
Variable  u  is  real,  while  A  can  be  complex.  For  the  non-central  chi-squared  RV  s  defined  in 
equation  (4),  the  mixture  function  is  given  by  equations  (23)  and  (9)  as 


U 

cs  (u,  A)  =  J  dx  ps  (x)  exp  (A  x) 


2  a  m 


2 


f  dx  x(M  1)/2  exp 

f 

-  X 

\-x) 

\ 

I M- 1 

(  r\ 

m-qx 

2 

J 

0 

V 

/2cr  K 

) 

l  J 

for  u  >  0, 


(24) 


where 


2  a  a 


(25) 


Upon  making  the  substitution  x  =  t2  in  equation  (24),  there  follows 


cs(u,A)  =  c(u,A;K,m,a)  =  rJ^  f  dt  tM  exp 

/r  m  J 


h 


2  M 


exp 


i-A- 

v  h  jj 


0 

r  h 


—  A 


f  mt^ 


r  n  i— 

qM[haU 


v2cr  jj 

1 


1  M-\ 


W  J 


for  all  A* - -,  a  >  0, 

2cr2 


(26) 


which  follows  directly  from  definitions  (1)  and  (20),  and  auxiliary  definitions 


AT  m  1 
M  -  — ,  r  -  —  =  — 
2  <7  cr 


2  <u(k)2  ,  h  =  h(A,<j)  =  Vl-2<J2A. 

V  k=\ 


(27) 


The  special  case  of  A  = - -  is  presented  later;  see  equation  (57). 

2a 

As  partial  checks  on  relation  (26),  there  follows  cs  (0,2.)  =  0 ,  upon  use  of  qM  {a, 0)  =  0. 
Also,  cs  (u, 0)  =  qM(r,  Vm /cr)  =  1  -  es (u),  which  is  in  agreement  with  equation  (20).  Finally, 
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Cs(+°M)  reduces  to  the  MGF  in  equation  (7),  using  qM(a,+oo)  =  1.  Result  (26)  has  been 
checked  numerically  against  the  integral  relation  (24). 

Another  quantity  of  interest  is 

8  u 

cs  (u,  A,l)  =  —  cs  ( u ,  A)  =  J  dx  x  ps  (x)  exp(A  x),  ^ 

by  reference  to  equation  (24).  By  making  the  change  of  variable  x  =  t2  and  referring  to  the  top 
line  of  equation  (26),  there  follows 

c.(M,D=-!2Ld7)  A  1  ,Y|, 

a2  mM~l  J  P  1  7m-i  ~T  ■  (2S 

o  \  v zcr  \<y  J 


At  this  point,  the  integral  result 


\dxxM+2  exp  -  £JL.  (ax)  =  4 


2M+4  )exp  -  j-  I  qM\  —,bp 

P  \2P  J  VP 


bM  f 

—  exp  - 
P  l 


p2b2 


a  I M  (a b)  +  b  p 2  I (ah)] 


is  used,  with  the  end  result  (after  some  manipulations) 


cs(u,A,\)  =  —  exp  - 


r2)\r2+2Mh2  f  r2  )  (r 

T  — US — exPkrr 


rM- 1  exP  2  j\-r  Im(vv)  + h2  v  IM](rv}]\  iorA^—j, 


j,  u  >  0,  m  >  0, 


K  Ti%  i  X  I — 

M  =  T’  r  =  ~  =  ~\I<V(k)2  ,  h  =  h(A,<T)  =  Jl-2a2A,  v  =  ^~. 


The  result  in  equation  (31)  has  been  checked  numerically  against  the  integral  relation  (29). 
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SADDLEPOINT  LOCATION  FOR  MGF  (7) 

The  logarithm  of  MGF  //s  (A)  in  equation  (7)  is 

Is W  =  log  As W  =  ~  lo80 -2a2A)  +  -m  4t  for Re(A)  < 

z  1  —  2(7/1  zcr 

The  derivative  of  this  quantity  is 

z's(A)  =  Ka2  y  +  m2y2,  y  =  1 

The  saddlepoint  (SP)  location  T0  is  given  by  the  solution  of  the  equation 
x  =  l's(A0)  =  Kcr2y0  +  m2y20. 


(33) 


(34) 


(35) 


where  x  is  the  field  point  of  interest  in  the  PDF  ps.  Solving  the  quadratic  equation  (35)  fory„ 
and  then  solving  the  second  term  of  equation  (34)  for  the  SP  gives 


K=K(x)  = 


la2 


1- 


2m2 


K2<j 4  +4  m2x-Ka2 


2cj2 


+  4m  x  +  Ka 


2x 


.(36) 


This  is  an  explicit  relation  for  the  SP  location  in  terms  of  field  point  x.  As  x  tends  to  zero,  A0 
tends  to  -oo,  while  as x  gets  large,  Aa  tends  to  l/(2cr2)-.  The  actual  saddlepoint  approximation 
to  the  non-central  chi-squared  PDF  at  field  point  x  is 


Ps(x) 


As(^0)exp(-T0x) 

pxx'M.) 


K o  = 


(37) 


where 


Z"(A)  =  2KcrAy2  +4m2a2y3 


(38) 


from  equation  (34). 


7 


EVALUATION  OF  FUNCTION  QM(a,b)  FORM  A  SMALL  HALF-INTEGER 


The  evaluation  of  QM(a,b),  defined  in  equation  (1),  is  simpler  when  M  is  a  half-integer, 
that  is,  when  the  number  K  of  RVs  added,  is  odd;  see  equations  (4)  and  (10).  For  M=  1/2,  the 
Bessel  function  in  equation  (1)  is  given  by  reference  8,  equation  10.2.14,  as 


7-i /20)  = 


f  n  A 


\JCZ  j 


1/2 


cosh(z). 


(39) 


Substitution  in  equation  (1)  immediately  leads  to  the  result 


01/2  (a,  b)  =  0(a  -b)  +  0(-a  -  b ),  gU2  (a, b)  =  1  -  Qin  (a, b)  =  0(6  -  a)  -  <J>(-tf  -  6),  (40) 


where 


O(z)  =  (2/r)'1/2  ^dt  exp(-f2/2)  =  (2tt)  1/2  jr#  exp(-tz/2)  (41) 

-co  — z 


is  the  normalized  Gaussian  CDF.  The  values  of  Qj[a,b)  for  M=  1 .5,  2.5,  3.5  ...  can  then  be 
found  from  reference  4,  equation  (3),  top  line,  which  does  not  require  M  to  be  an  integer.  The 
required  values  of  the  Bessel  functions  are  available  in  reference  8,  equation  10.2.13.  The  end 
results  are 


Q3/2(a,b)  =  Q,/2(a,b)  +  J—  exp| 

n 


f  a2+b2^ 


sinh(a  b ) 


V  A 

f  J  ,  L2\ 


a 


Qsn(a’b)  =  Qm(.a,b)  +  exp) 

n 


a1 +b2 


a  b  cosh(a  b)  -  sinh(a  b) 


V  z 
{  _2  ,  .2  \ 


a 


(42) 


Q7l2(a,b)  =  Q5/2(a,b)  +  J—  exp 

n 


a 2  +b 2 


(3  +  a 2  b 2 )  sinh(er  b) -3  a  b  cosh(a  b) 


a ' 


The  quantity  qM(a,b)  appears  in  equations  (26)  and  (31).  Then,  numerically,  the  second 
result  in  equation  (40)  is  more  appropriate,  and  the  top  line  of  equation  (42)  can  be  modified  to 


2 

{  a2+b2\ 

-  exp 

— 

n 

V  2  J 

sinh(a&) 


a 


(43) 


along  with  obvious  changes  to  the  other  two  relations.  This  procedure  avoids  the  subtraction  of 
Qm  from  unity  and  retains  more  accurate  numerical  results. 
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All  three  results  in  equation  (42)  have  removable  singularities  at  a  =  0.  The  special  values 
there  are 


Q3i2(0,l>)  =  2Q(-b)  +  J—  exp 


f  b ^ 


n 


\  2  ) 


b, 


Qs/2 (0,b)  =  03/2 (0.*)  +  -.  -  exp 


[ 2 

{  b2) 

J—  exp 

— 

V  n 

<  2  J 

[2 

f  b2) 

t  —  exp 

— 

V  n 

^  2  J 

(44) 


3  1 

15 


However,  the  last  two  results  in  equation  (42)  have  numerical  problems  when  product  a  b  is 
near  zero,  due  to  the  differences  of  nearly-equal  quantities  in  the  numerators.  Power  series 
expansions  are  required  for  the  last  two  results  in  equation  (42);  they  are 


z3  z5  z7  z9  z"  z13 

z  cosh(z)  -  sinh(z)  =  —  +  —  H - h - + - h - . 

3  30  840  45360  3991680  518918400 


z5  7. 7  7 

(3  +  z2 )  sinh(z)  -  3  z  cosh(z)  s  — -  + - + 


.n 
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.15 


(45) 


■  +  - 


15  210  7560  498960  51891840  7783776000 


A  MATLAB  program,  QM_half_integer(M,a,b),  has  been  written  that  incorporates  all  these 
features,  at  least  up  to  M-  7/2.  Although  a  recursion  based  on  reference  4,  equation  (3),  may 
seem  attractive,  it  is  not  useful  for  small  a,  as  may  be  seen  from  the  low-order  results  in  equation 
(42)  above,  where  successively  higher  powers  of  parameter  a  appear  in  the  denominator.  The 
recursion  inevitably  encounters  differences  of  like  quantities,  which  are  then  divided  by  very 
small  quantities,  resulting  in  severe  losses  of  significance  in  the  final  results.  Nothing  in  a 
recursion  can  avoid  the  required  switch  to  a  power  series  (45)  for  small  arguments. 


SERIES  EVALUATION  OF  qA1(a,b)  FORM  INTEGER 
The  series  expansion 


00  A *  k+M-l  D H  2  l2 

e„ («■  6)  =  expM  -  B)  £  A  =  =  (46) 

k=o  k.  „=0  n\  z  z 

was  derived  in  reference  4,  equation  (8),  by  expanding  the  Bessel  function  I M  , (ax)  in  equation 
(1)  in  a  power  series  and  integrating  term  by  term.  By  expressing  the  inner  sum  on  n  as  the 
difference  between  a  sum  to  infinity  and  a  sum  from  k  +  M  to  infinity,  equation  (46)  easily  yields 
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(47) 


The  last  form  involves  only  one  infinite  sum,  all  the  terms  of  which  are  very  amenable  to 
evaluation  by  way  of  recursions.  Also,  for  a  and  b  real,  every  quantity  is  nonnegative. 

However,  if  a  or  b  is  complex,  positive  and  negative  quantities  will  occur  in  equation  (48),  and 
all  significance  can  be  lost  in  the  final  sum.  Equation  (27)  indicates  that  if  A  is  complex,  then 
h  =  h(A,cr)  will  be  complex,  and  the  gM(a,b )  term  in  equation  (26)  will  have  to  be  evaluated 
for  complex  arguments  a  and  b. 

When  |/?|  is  large  relative  to  \a\,  a  useful  alternative  to  the  above  series  expansions  is 
afforded  by  equation  (4)  of  reference  4: 

QM(a,b)=  -  exp  -  —  (b-a)2  -  exp (-ab)  I„+l_M(ab).  (49) 

\aj  V  2  Jn=o\bJ 


The  magnitudes  of  both  terms  in  the  sum  eventually  decrease  with  n;  see  reference  8,  page  428. 
The  useful  combination  in  equation  (12)  is  encountered  again  in  this  form. 

Conversely,  when  |o|  is  larger  than  \b\,  equation  (5)  of  reference  4  yields 

f  h^\M  (  1  \  »  ( u\n 

qM(a,b)=  -  exp  --(b-a)2  £  -  exp(-a&)  I„+M(ab).  (50) 

\a)  V  2  J  „=0  \a) 
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JOINT  PROBABILITY  DENSITY  FUNCTION  OF  J- 1  LARGEST  RANDOM 
VARIABLES  AND  THE  SUM  OF  THE  REMAINDER  OF  A 
SET  OF  N  INDEPENDENT  RANDOM  VARIABLES 


The  joint  PDF  of  the  J-  1  largest  RVs  and  the  sum  of  the  remainder  of  a  set  of  N 
independent  RVs  {sn}  is  presented  in  reference  7 — in  particular,  see  equations  (15),  (28),  (38), 
and  (47),  for  J=  2,  3, 4,  and  5,  respectively.  The  PDFs  {p„(x)}  of  the  N  independent  RVs  {s„ } 
are  arbitrary  in  all  these  cases. 


NON-CENTRAL  CHI-SQUARED  RANDOM  VARIABLES 

Let  the  «-th  independent  RV  s„  be  given  by 


s„=Zg»(*)2  f°r  n  =  l:  N, 


(51) 


k= 1 


where  independent  Gaussian  RVs 

g„ (k)  =  Normal{//„ (*), cr„ }  for  k  =  l:Kn,  n  =  l:N. 


(52) 


Notice  that  standard  deviation  an  is  independent  of  k.  Then,  as  shown  in  equation  (9),  the  PDF 
ofRV  s„  is 


(u)  = 


2  a 


\r*j 


exp 


^  t2  +r2  ^ 

n  n 


IMn-\(rn  O  for  U  >  °»  mn  >  0>  (53) 


where 


,,  Kn  (K 

Mn  mn  = 


\1/2  /— 

X  Pn  (kf  »  rn=~’  tn=—  for  «  =  1  I  N. 


V*= i 


(54) 


Also,  from  equation  (26),  the  mixture  function  of  RV  s„  is 


cn(u,A)  = 


where 


1 


hn(A)2M» 


exp 


KW 


J) 


Vn  KW 


hMY 


for  all  X  #  — X.,  u  >  o,  (55) 
2<j„ 


hn(X)  =  h(X,an)  =  yjl-2a2  X  for  n  =  1 :  N. 


(56) 
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For  A  =  /?n  (A)  =  0  and  equation  (55)  reduces  to 

2cr 


r  1 1 

(  r\ 

VM 

c„ 

=  exp 

rn  -  , 

l  2(J«; 

v  2  J 

U  ^ ; 

l  aJ 

,  u  >  0,  mn>  0. 


(57) 


For  mn  =  0,  this  reduces  further  to 


1 


V  lGn  ) 


_1_ 

M\ 


(  \ 


M 


V2^7 


,  u  >  0. 


(58) 


Interest  is  centered  here  on  the  situation  where  Na  of  the  K a  sets  of  N  independent 
Gaussian  RVs  in  equation  (52)  have  a  common  mean  and  standard  deviation,  namely, 

g„(k)  =  Normal^ (k),aa}  for  k  =  \:Ka,  n  =  \:Na, 


(59) 


while  the  remaining  Kb  sets  of  Nb  =  N  -Na  independent  Gaussian  RVs  have  a  different 
common  mean  and  standard  deviation: 


g„(&)  =  Normal{//6 (k),ah)  for  k  =  \:Kb,  n  =  Na+l:N. 
Then,  from  equation  (8), 


m„  = 


( K*  V/2 

l>a(*)2 

V*= 1 


>  = 


(Kb 

Z  Mtikf 

V*=i 


(60) 


(61) 


The  corresponding  PDF  and  mixture  function  of  RVs  {s„},  n  =  1 :  Na,  are,  respectively, 


<  ^  ' 


Pa(u)  =  p(u;Ka,ma,cra)  = 


ca  (u,  A)  =  c(u.  A;  Ka,ma,cra ) 


2cr 


M- 1 


«  r°  ) 


exp 


2  2 


2a 


1 


ha(A)2M° 


exp 


f  rU 


1  — 


2  ha(AY 


'  Ma~  1 

r 


( 

\°°  J 


r  Ma 


JJ 


r°  h°{X)47? 


(62) 


for  u  >  0  and  n  =  1 :  Na ,  where 
K 


rB=- S  ha(A)  =  h(A,aa)  =  yll-2a2A. 
2  C5" 


(63) 


For  RVs  {s  n},n-Na  + 1 :7V,  replace  all  the  subscripts  a  by  b  in  equations  (62)  and  (63). 
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JOINT  PDF  OF  LARGEST  RANDOM  VARIABLE  AND 
SUM  OF  REMAINDER 

The  first  case  of  interest  is  given  by  equations  (15)  and  (16)  of  reference  7;  namely,  the  joint 
PDF  is,  for  J=  2, 


?2  0,  ,z2)  =  7T  |  dX  exp(-/l  z2 )  P(z, ,  A)  £  , 

l2”  c  tt  c„(zx,A) 


(64) 


where  product 

P(zj  ,A)  =  Y[c„  (z,  ,A)  =  ca  (z, ,  A) N°  cb  (z, ,  A)Nb  (65) 

»= i 


and  sum 

y  P.Ah)  _N  PaiZx)  LV  pfc(^) 

»= 1  C„(Z15^)  °  ca(zx,A)  b  cb(z1,A)' 


(66) 


Bromwich  contour  C  in  equation  (64)  runs  from  -  zoo  to  +  zoo  in  the  complex  /l  plane.  For 
numerical  accuracy,  this  contour  C  can  be  moved  so  as  to  pass  through  the  (unique  real)  SP  As 
of  the  integrand  of  equation  (64);  of  course,  this  SP  location  depends  on  the  particular  two- 
dimensional  field  point  (zx  ,z2)  of  interest,  that  is,  As  =  As  ( z , ,  z2 ).  This  SP  location  must  be 
recalculated  every  time  the  two-dimensional  field  point  is  changed. 


JOINT  PDF  OF  TWO  LARGEST  RANDOM  VARIABLES  AND 
SUM  OF  REMAINDER 

This  result  is  available  from  reference  7,  equations  (28)  through  (30);  namely,  the  joint  PDF 
is,  for  J=  3, 


q3(z  1  ,z2,z3)  =  U(Zl  -z2) 


1 

i2n 


J dA  exp (-Az3)  P(z2,A )  (Sj 

c 


S2 


•U 


(67) 


where  product 

P(Z2,A)  =  Y[  C»(Z2’A)  =  Ca(Z2’A)N°  Cb(Z2,A)Ni 


N 


«= 1 


(68) 


and  sums 
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(69) 


Si  =  NaMhL  +  NiJ^L'  s  N  E^hL+Nb  Edhl  , 

ca(z2,X)  ch(z2,A)  “  cfl(z2,A)  cb(z2,A) 

5  =iV  pAeApJeA  i  jy  Pb(zi)pb(z2 ) 

3  *  Ca(z2,A)2  6  CA(z2,yl)2 

The  comments  under  equation  (66)  are  again  directly  relevant,  except  that  the  field  point 
(zi,z2,z3)  in  equation  (67)  is  now  three-dimensional. 

The  extensions  to  larger  values  of  J  are  given  in  reference  7,  equations  (38)  through  (42)  for 
J  =  4,  equations  (47)  through  (50)  for  J  =  5,  and  equation  (51)  for  the  general  case. 
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EVALUATION  OF  THE  Q  FUNCTION  FOR  ALL  REAL  ARGUMENTS 


In  this  section,  attention  is  limited  to  the  case  M  =  1.  Then,  equation  (1)  reduces  to 


x2  +  a2  ^ 


I0(ax). 


(70) 


Although  this  function  is  analytic  in  each  of  the  variables  a  and  b,  initial  interest  here  is  in 
numerical  evaluation  of  Q(a,b)  only  for  nonnegative  real  a,b  for  all  ranges  of  their  values.  One 
of  the  problems  associated  with  this  evaluation  is  that  underflow  and/or  overflow  frequently 
occur  for  large  arguments,  and  all  significance  is  lost.  A  method  will  be  presented  that  isolates 
the  underflow/overflow  problem  and  allows  for  the  calculation  of  Q  or  log(0  for  all  arguments. 
Even  when  Q  underflows  or  overflows  the  capability  of  the  computer,  log(0  still  retains 
significance.  This  situation  is  analogous  to  the  gamma  function  T(x),  where  both  F(x)  and 
log(T(x))  are  furnished  by  a  computer,  to  circumvent  overflow. 

In  addition,  the  complementary  Q  function  now  becomes 


x2  +a2^ 


I0(ax). 


(71) 


When  Q  approaches  1  so  closely  that  all  significance  is  lost,  the  complementary  q  function  still 
retains  significance.  Again,  since  q  could  also  underflow,  there  is  interest  in  computing  both  q 
and  log(g)  for  all  a,  b.  This  situation  is  similar  to  that  where  both  the  error  function,  erf(x),  and 
the  complementary  error  function,  erfc(x)  =  l-erf(x),  are  furnished  by  a  computer.  The  Q 
function  is  an  exceedance  probability,  while  the  q  function  is  a  cumulative  probability  of  a 
normalized  Rician  random  variate,  as  noted  in  an  earlier  section. 


In  reference  3,  equation  (74),  the  following  integral  result  was  presented: 


1 

[  dx  —  exp 

'  n 


a 


1  -  P  cos(x) 


f  1  _  1  _ ^ 

=  2  Q  —  yja  (1  -  w),—  yla  (1  +  w) 
\u  u  J 


-expf 


a 


\-H2 


a  (3 

W?* 


(72) 


for  ft  <  1,  with  u  =  yjl-  jB2 .  By  making  appropriate  changes  of  variables,  this  result  can  be 
manipulated  into  the  following  forms  (reference  5,  equation  18) 


Q(a,b)  =  I  +  E,  q(a,b)  =  1-  Q(a,b)  for  a<b, 


(73) 


and 
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(74) 


q(a,b)  -  I  -  E,  Q(a,b)  =  1  -  q(a,b)  for  a>b. 


where 


I  =  — —  f  dx  exp - — -  c  r  ( a  b  )  _  lab 

ln  0  l  l-c2cos(x)J’  1  2  (o2  +  b1)  ’  2  a2  +  b2  ’ 

P  1  f  a2+Z>2^| 

£  =  -exp - —  /o(ai). 


A  useful  alternative  to  the  latter  term  is 

p  f  (a-b)2)  1 
E  =  exp - 2  exP(~*  b )  4  («  *), 


in  thatthe  product  exp(-x)  70(x)  is  very  well  behaved  for  all  real  nonnegative x,  decaying  as 
1  -<J27rx  for  large  x.  In  fact,  this  product  is  denoted  as  besseli(0,x,l)  in  MATLAB.  Thus,  all 
the  underflow  problems  associated  with  E  are  isolated  in  the  leading  exponential. 

At  the  same  time,  integral  /  in  equation  (75)  can  be  expressed  as 


/  =  exp| - — —  — L.  f  dx  ext/-  c‘  Cz  1  +  cos(x)  ) 

l  1  +  c2  J  2n  *  P[  l  +  c2  l-c2cosfx)l’ 


cos(x) 


for  which  the  maximum  value  of  the  integrand  is  1  at  *  =  *  Now,  by  employing  the  definitions 
in  equation  (75),  this  expression  takes  the  form 

/  =  expf-  J_  ]dx  expf_  £ )2  J+cosQQ  \ 

v  ^  J2^0  (_  a2  +b2  l-c2cos(x)J 

Again,  the  exact  same  leading  exponential  accounts  for  all  the  underflow  problems  associated 

with  /.  The  integrand  of  equation  (78)  increases  monotonically  over  the  interval  (0 ,tt),  ending 
up  at  value  1 ,  and  cannot  lead  to  underflow. 

For  computational  purposes,  define  the  parameters 

(79 


Then,  from  equations  (78)  and  (76),  there  follows 
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(80) 


which  will  not  underflow  for  any  a,b.  The  complementary  quantities  in  equations  (73)  and  (74) 
are  then  available  upon  subtraction  of  equation  (81)  from  1.  At  least  one  of  the  four  quantities 

Q,  q,  log(0,  log(gO  (83) 

will  retain  significance  for  all  values  of  a,b. 

A  MATLAB  routine  was  written  for  evaluation  of  the  four  quantities  in  equation  (83),  and  is 
called  according  to 

[Q,q,Q  log,  q  log]  =  Qq  log  (a,  b )  (84) 

The  time-consuming  portion  of  this  procedure  is  the  numerical  evaluation  of  the  integral  in 
equation  (82).  The  major  problem  is  that  the  monotonic  integrand  in  equation  (82),  namely, 


g(x) s exp 


^4 


1  +  cos(x) 
l-tj  cos(x) 


(85) 


varies  from  g(0)  =  exp(-  2  t4  /(I  - 13 ))  =  exp(-2  a  b )  to  g(;r)  =  1,  and  can  do  so  in  a  very  narrow 
interval  in  x.  This  can  sometimes  cause  the  routine  in  equation  (84)  to  expend  a  great  deal  of 
time  achieving  a  stable  estimate  of  the  integral  required  in  equation  (82).  Integrand  g(x)  does 
not  have  a  peak  anywhere  in  (0,tz);  however,  g(x)  can  have  a  narrow  transition  range  for  its 
values,  which  impedes  quick  evaluation  of  the  integral.  This  sampling  problem  has  been 
circumvented  by  resorting  to  MATLAB’s  quad  or  quadf  routines,  which  employ  an  adaptive 
Simpson’s  rule  that  resorts  to  finer  sampling  where  the  integrand  varies  most  rapidly. 
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In  an  effort  to  ameliorate  this  rapid  transition  of  g(x),  an  integrand  of  similar  shape  was 
constructed,  namely, 


h(x)  =  exp 


-tA 


1  +  cos(x) 
1  -  cos(x) 


=  exp(/4 )  exp 


2  /, 


y  1  -  cos(x) 


(86) 


Its  integral  is  available  in  closed  form  (reference  3,  equation  (75)): 


|  dx  h(x)  =  -j—  exp(t4 )  J  exp 


2  t. 


1  -  cos(x) 


=  ]-  exp(t4 )  erfc(^7 )  =  ^  erfc ),  (87) 


where  the  last  function  is  an  available  MATLAB  routine.  The  function  erfcx(x)  is  well  behaved 
for  all  positive  real  x  and  decays  as  1 / (V/r  x)  for  large  x. 

The  pertinent  integral  in  equations  (80)  through  (82)  can  now  be  expressed  as 


1 

2  n 


71 

|  dx  g(x) 
0 


~  f  dx  [g(x)  -  h(x)\  +  ^  erfcx(>/J7), 
2  n  0J  2 


(88) 


leaving  the  major  numerical  problem  as 


( 


w(x)  =  g(x)  -  h(x)  =  exp 


l  +  cos(x)  ^ 
1  -/3  cos(x)y 


f 

-exp 

v 


U 


1  +  cos(x)^ 
1  -cos(x). 


(89) 


The  difference  function  w(x)  has  the  special  values 
w(0)  =  exp(-2  a  b),  w(7r  /  2)  =  0,  w(x)  =  0, 


(90) 


indicating  that  a  major  portion  of  function  g(x)  has  been  accounted  for,  by  means  of  function  h(x) 
defined  in  equation  (86).  Integral  (89)  is  also  a  good  candidate  for  MATLAB ’s  quad  or  quad t 
routines. 


EVALUATION  OF  Q{a,b)  FOR  COMPLEX  ARGUMENTS 


The  integrand  in  equation  (85)  can  be  extended  to  the  complex  z-plane  by  analytic 
continuation: 


( 


g(z)  = exp 


1  +  cos(z) 

1  - 13  cos(z)  y 


(91) 


18 


The  function  g(z )  has  an  essential  singularity  (ES)  where 


cos(z0)  =  l/f3,  z 0  —  acos(l/?3 ),  t3 


lab 
a2  +b2 


(92) 


However,  this  complex  value  z0  is  only  the  principal  value  location  of  the  ESs  of  g(z).  The 
totality  of  ESs  of  g(z)  are  located  at 

z0+2nn  and  -z0+l7rn,  n  =  0,  ±1,  ±2,....  (93) 


Principal  value  zG  has  a  real  part  that  always  lies  in  the  interval  [0,;r].  Therefore,  the  only  ESs 
that  can  affect  the  integration  over  [0,zr]  in  equation  (82)  are  those  at 

±  zc  and  ±  zQ  +  2  n.  (94) 


If  input  arguments  a  and  b  are  real,  positive,  and  unequal,  then  t3  is  real  and  less  than  1,  and 

z0  =  acos(l/ >3 )  =  /  acosh(l/ 13 ) ,  (95) 

which  lies  on  the  imaginary  axis  of  the  z-plane.  Another  ES  lies  at  -  z0 .  Therefore,  the  path  of 
integration  in  equation  (82)  starts  out  perpendicular  to  the  line  between  the  two  ESs  on  the 
imaginary  axis.  This  path  is  the  best  way  of  avoiding  the  neighborhoods  of  the  two  closest  ESs, 
and  yet  begin  the  integration  at  z  =  0.  The  two  other  ESs  in  equation  (94)  are  located  above  and 
below  2n  in  the  z-plane  and  do  not  significantly  affect  the  integrand  in  the  interval  [0 ,zr]. 

On  the  other  hand,  when  either  input  argument  a  and  b  is  complex,  principal  value  location 
z0  in  equation  (92)  can  have  a  real  part  that  approaches  n.  Then,  the  real  part  of  2  n-zQ  also 
approaches  n,  thereby  possibly  leading  to  two  ESs  close  to  the  end  of  the  path  of  integration  at 
z  =  it.  To  try  to  stay  away  from  the  immediate  neighborhoods  of  these  ESs,  the  following 
modification  of  the  path  of  integration  in  the  z-plane  is  suggested.  First,  draw  a  line  Lx  between 
z0  and  -  z0 .  Then,  start  the  integration  from  z  =  0  along  a  line  Px  perpendicular  to  Lx .  Second, 
draw  a  line  L2  between  z0  and  2 n-za.  Finish  the  integration  at  z=  n along  a  line  P2 
perpendicular  to  L2 .  This  approach  attempts  to  avoid  the  neighborhoods  of  the  ESs,  wherever 
they  may  be  located.  The  lengths  of  the  two  perpendicular  lines  Px  and  P2  could  be  taken  as 
7if2\  this  guarantees  that  the  projections  of  the  two  perpendicular  lines  Px  and  P2  on  the  x-axis 
will  not  extend  beyond  each  other,  yet  allows  the  path  of  integration  to  get  a  sufficient  distance 
away  from  the  ESs.  Finally,  join  the  ends  of  the  two  perpendicular  lines  Px  and  P2  with  another 
line,  thereby  completing  the  modified  contour  between  z  =  0  and  z  =  n. 
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This  procedure  has  been  found  to  be  adequate  for  some  complex  pairs  of  arguments  a,b,  but 
not  for  all  arguments.  When  the  magnitudes  of  arguments  a,  b  get  large,  some  pairs  of  arguments 
cause  large  variations  in  the  exponential  functions  involved,  and  all  significance  can  be  lost. 

ALTERNATIVE  EVALUATION  OF  THE  Qu(a,b)  FUNCTION 

The  QiJ^a.b)  function  can  be  expressed  in  terms  of  the  Q(a,b )  function  according  to 
(reference  4,  equation  (3)): 

(  (n f  h\k 

Qm  (a,  b )  =  Q(a,  b)  +  exp - — —  £  -  exp(-a  b )  I k  ( a  b).  (96) 

V  2  )  k= i  W 

Again,  the  same  underflow  factor  as  in  equations  (76)  through  (81)  appears  in  the  additive  term; 
additionally,  the  well-behaved  combination  exp(-x)  1  k  (x)  is  encountered.  Therefore,  a  simple 
modification  to  equations  (81)  and  (82)  enables  evaluation  of  QM(a,b )  via  equation  (96).  The 
apparent  danger  of  small  a  in  the  b/a  term  is  partially  compensated  for,  by  the  fact  that 

—  as x  — »  0.  (97) 

x*  2 k  k\ 

However,  for  complex  a,b  if  the  sequence  (exp(-a&)  Ik  ( ab )}  has  been  determined  by  means  of 
a  recursion  using  subtraction,  form  (96)  will  lose  significance  for  small  a*  0.  Then,  it  will  be 
necessary  to  employ  a  power  series  expansion  of 

M-\  f  f.\k 

£  -  w 

m 

about  a  =  0.  The  case  of  a  =  0  must  be  taken  care  of  separately. 

A  routine  for  the  evaluation  of  the  four  quantities 

QM(a’hX  9a/ (A 6)  =  1“ 2m («,£)>  log (QM(a,b)),  log (qM(a,b))  (99) 

was  written  and  made  available  with  the  call 

[Q,  q,  Q  log,  q  log]  =  Qq\ogM(M,a,b).  (100) 

At  least  one  of  the  four  quantities  in  equation  (99)  retains  significance  for  all  arguments.  This 
routine  can  be  extended  to  complex  arguments  a, b,  but  is  not  trustworthy  for  all  complex 
arguments. 
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EVALUATION  OF  THE  0„(a,b)  FUNCTION  FOR  COMPLEX  ARGUMENTS 
BY  MEANS  OF  INTEGRATION  ON  THE  REAL  AXIS 


The  QM(a,b )  function  was  defined  in  equation  (1)  as 


uu 

Qu(a,b)=  \dxx 


exp - - —  /«.,(«). 


(101) 


The  integrand  is  analytic  for  all  finite  complex  a,  thereby  leading  to  the  QM  (a,  b)  function  being 
analytic  in  both  variables  a,b  for  all  finite  complex  values.  The  complementary  QM  function 
was  defined  as 


2  2 

C  X  X*“  |  Cl*" 

qM(a,b)  =  l-QM(a,b)=  f  dx  x  -  exp - - -  1  M_ ,(«). 

i  \  a  2 


(102) 


By  expanding  the  Bessel  function  in  equation  (101)  in  a  power  series  and  integrating  term  by 
term,  an  alternative  expression  is  obtained  (see  reference  4,  equation  (8)): 


^  ,  l\  a2+b2' ^(a2  2)k  k+^~l  (b2  2)” 

QM(a,b)  =  cx p - —  X ~ 7j  ~  L  — , — 

l  2  i  to  k\  n  n\ 


(103) 


This  form  allows  for  a  useful  observation,  namely,  that  the  QM  ( a,b )  function  is  even  in  both 
arguments.  That  is, 

Qm  («> b )  =  Qm  (~a> b )  =  Qm  (a~b)  =  Qm  i-o-b),  ( 1 04) 

even  when  arguments  a  and  b  are  complex.  Advantage  will  be  taken  of  this  fact  later.  However, 
direct  use  of  expansion  (103)  for  numerical  evaluations  leads  to  underflow  and  overflow  for 
moderately  large  values  of  arguments  a  and  b,  and  to  cancellation  and  loss  of  significance  for 
complex  a, b.  Accordingly,  an  alternative  representation  of  the  QM  ( a,b )  function  is  sought  that 
will  allow  its  accurate  evaluation  for  large  complex  arguments  a, b. 

By  repeated  integration  by  parts  on  equation  (101)  (see  reference  4,  equation  (4)),  another 
representation  is  achieved: 

OO  2  1  2 

QM(a,b)  =  exp(-A- B)  V  tk  Ik{u),  where  A  = — ,  B  =  — ,t  =  —,u~ab.  (105) 

k=\-M  2  2  b 

This  expression  can  be  developed  as  follows: 
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(106) 


G„(o,6)  =  expM-B)/“  jr  (**'";il(»)  =  expM-B)(,-"X/'7ytl^(U) 

fc=l-A/  y=0 

1  71  oo 

=  exp (-A  -  B)  tl~M  —  |  dx  exp[a  cos(x)]  ^  t1  cos[(y  + 1  -  M)x\ 


The  sum  on  j  can  developed  into  a  more  useful  form  as 


-  X  fi  {exPl>‘  U  + 1  “  M)  x]  +  exp[-/  (J  + 1  -  M)  x] } 

2  ;=o 

1  00  1  00 

=  —  exp[-i  x  (M  - 1)]  ^  [/  exp(z  x)]7  H —  exp[z  x  ( M  - 1)]  ^  [/  exp(-z  x)]7 
2  y=o  2  j= o 

=  I + 1  exppx(Af-l)]  &r  .  ,  <  (hat is> y  < |A| 

2  1  -  /  exp(z  x)  2  1-/  exp(-i  x) 


(107) 


Substitution  in  equation  (106)  yields 


~  ,  ,N  ,  A  Dsx-m  1  r  ,  r  ,  vJ  exp[zx(A/-l)]  exp[-zx(M-l)] 

Qm  b )  =  exp(-^  -  B)  t  —  J  dx  exp[zz  cos(x)]  - : - — — 

2/r  -  \^l-/exp(-zx)  l-/exp(zx) 

=  expM -  B) I  f  A exp[w  cos(x)]  60s[(M-l)xl-l  cos(A7x) 
n  *  l-2/cos(x)  +  / 


(a-iO’U&V  1  L  r  ,  w,  '  ' 

=  -  exp - —  —  ax  exp[-ao(l-cos(x))J - - — 

l  2  la)  2^i  1  lab 


cos  (M  x) cos((M  - 1)  x) 


(108) 


-  -  +  -  -cos(x) 
2\b  a 


for  ja|  <  |/>|,  where  the  definitions  in  equation  (105)  have  been  utilized. 

In  an  exactly  analogous  fashion,  the  expansion  for  complementary  function  qM(a,b )  in 
reference  4,  equation  (5),  can  be  developed  into  an  integral  representation.  The  end  result  is 
identical  to  that  in  equation  (108),  except  for  the  minus  sign  out  front  and  the  alternative 
restriction  now  that  |#|  <  |a|.  Therefore,  the  situation  can  be  summarized  as  follows:  define 
integral 


I  =  exp  - 


a  )  2  7t 


n  cos (M  x)  —  cos((M  - 1)  x) 

\dx  exp[-a6(l  -  cos(x))] - -  (109) 

J  1  a  b  i 

0  1  -  +  -  —  cos(x) 

21  b  a  J 


Then  (see  also,  reference  9,  equations  (C-26)  and  (C-27)) 


22 


(110) 


qM(a,b)  =  +I  if  |Z>|<|a|, 
QM(ci,b)  =  -I  if  |a|<|4 


The  main  issue,  then,  is  the  evaluation  of  the  integral  I  in  equation  (109)  for  complex  a,b. 
Extend  the  integrand  of  equation  (109)  into  the  complex  z-plane  by  analytic  continuation  and 
define  numerator 


N(z)  =  exp[-<2&(l  -  cos(z))] 


|^cos(M  z ) - cos((M  - 1)  z)  j 


(111) 


This  function  is  entire  in  z;  that  is,  N(z )  has  no  singularities  anywhere  in  the  finite  z-plane.  Also, 
the  analytic  continuation  of  the  denominator  is 


Z>(z)  =  I 


(a  h  ' 
— +  — 
\b  aj 


-  cos(z). 


This  function  D(z)  has  simple  zeros  at  z  =  z0 ,  where 


cos(z0)  =  i 


a  b 
( b  a 


lab 

,  or  exp(iz0)  = 

t  b  a 


as  seen  from  the  first  line  of  equation  (108).  Therefore, 


iz0  =  i  2  n  n  ±  log 


r 


=  ilnn± 


log 


+  i  angle 


(112) 


(113) 


(114) 


where  n  is  an  arbitrary  integer,  positive  or  negative.  Here,  functions  log(z)  and  angle(z)  are  the 
principal  value  functions.  Finally,  the  locations  of  the  zeros  of  denominator  D(z)  are  at 


z„  =l7tn± 


Cu\ 


angle^— J  -  i  log 


(115) 


One  of  these  zero  locations  always  has  a  real  part  that  lies  in  [0 ,n],  which  is  the  region  of 
integration  for  I  in  equation  (109).  The  imaginary  part  of  this  zero  location  can  lie  above  or 
below  the  x-axis,  or  it  can  lie  right  on  the  x-axis  when  \a\  =  \b\. 

In  fact,  every  strip,  [nn,  (n+  l)zr],  of  the  z-plane  contains  a  zero  of  denominator  D(z). 
Furthermore,  some  pairs  of  these  zeros  can  lie  very  close  to  each  other.  For  example,  if  z0  is  near 
0,  then  the  zero  at  -  z0  is  also  near  0.  Or  if  ze  is  near  n,  then  2  n-z0  is  also  near  it.  Since  the 
zeros  of  D(z)  are  the  simple  poles  of  the  integrand  N(z)/D(z)  of  equation  (109),  direct  numerical 
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evaluation  of  /by  means  of  equation  (109)  will  run  into  difficulties  when  these  zeros  lie  close  to 
the  real  x-axis,  that  is,  j<j|  « |Z>|.  However,  equation  (109)  is  imminently  useable  when  \b/a\  is 
sufficiently  different  from  unity,  because  all  the  poles  of  N(z)/D(z)  are  sufficiently  removed 
from  the  neighborhood  of  the  path  of  integration  over  [0 ,7t\.  Furthermore,  equation  (109)  is  much 
more  attractive  than  equations  (101)  or  (102)  for  two  reasons:  no  Bessel  functions  need  to  be 
evaluated  and  the  range  of  integration  is  finite.  Also,  only  one  exp  function,  and  no  more  than 
three  cos  functions,  need  to  be  evaluated,  all  on  the  real  x-axis.  In  fact,  for  M  =  1 ,  only  one  cos 
function  is  involved;  for  M=  2,  only  two  cos  function  evaluations  are  required. 

At  first  sight,  it  might  appear  that  the  contour  of  integration  in  equation  (109)  could  be 
moved  away  from  the  pole  locations  of  N(z)/D(z),  by  moving  the  contour  upward  or 
downward,  into  the  half-plane  away  from  the  nearest  pole.  This  procedure  has  two  drawbacks. 
First,  the  cos  functions  in  equation  (109)  develop  exponential  behavior  because 

cos(x  +  i  y)  =  cos(x)  cosh(y)  -  i  sin(x)  sinh(y)  (116) 

for  nonzero  y,  thereby  leading  to  large  oscillating  behavior  of  the  integrand  with  x.  Second, 
when  the  poles  are  close  to  0  or  ?r,  it  is  not  possible  to  get  away  from  these  poles  because  the 
path  of  integration  must  start  and  stop  at  these  points. 

To  handle  the  case  where  jaj « |/[,  or  the  poles  are  close  together,  the  singularities  of  the 
integrand,  N(z)/D(z),  in  equation  (109)  will  be  subtracted  out  and  integrated  analytically.  To 
accomplish  this  task,  observe  from  equation  (113)  that 

a 
~b 

Therefore, 

h  M-2  (n2  —h2} 

cos  (Mza)  —  cos((M-l)zo)  = - — - .  (118) 

a  2b 


exp  (ikza)  = 


'a' 


\n  J 


k  (  u\k 


\aj 


;  cos  (kz0)  =  ^ 


+ 


( 


\a) 


(117) 


Reference  to  equations  (1 1 1),  (1 13),  and  (118)  then  yields 


N(za)  =  exp 


-  ab 


a  b  V 
— +  — 
\b  a)) 


a 


M-2 


(a2-b2) 


2b 


M 


-  exp 


-Aa-bf 


^  aM~2 


A 


(a2-b2) 


2b 


M 


•  019) 


Now,  express  the  integrand  of  equation  (109)  as 


The  integral  of  the  last  term  in  equation  (120)  is 


f  dx^^-  =  f  dx- 
.  D(x)  J  1 


D(x) 


N(z  )  2nab 

- =N(z0)  — — —  p  for  b*±a, 


a  b 
—  +  — 
b  a 


(121) 


-  cos(x) 


a 


where  polarity 

h1  if  H<I6I] 


p  = 


+  1  if  \b\  <  \a\ 


(122) 


Here,  the  integral  result  is  used 


n  1 

f  dx - 

*  a  -cc 


n 


cos(x)  4a  + 1  1 


if  ag[-l,l]. 


(123) 


where  both  square  roots  are  principal  value.  (The  branch  line  of  the  product  of  square  roots  is 
along  the  real  axis  of  the  complex  a  plane,  between  a  =  -1  and  a  =  1 .  The  denominator  of 
equation  (123)  cannot  be  replaced  by  principal- value  square  root  yja2  -1.) 


Combining  the  results  in  equations  (119)  and  (121)  yields 


]  dx  n  exPl 

0 


D(x) 


f  i 

y-6)J 


M- 1 


a 

\bj 


P- 


(124) 


Finally,  use  of  this  result  and  integrand  breakdown  (120)  yields,  for  the  second  component  of 
equation  (109),  the  very  simple  result 


,  1 
I2=~P- 


(125) 


The  first  component  of  equation  (109)  is  given  by 


/,  =  exp 


f 


0 a-b ) 


2 


\aj 


1  n 

-  f  dx 

In  J 


N(x)-N(z0) 

D(x) 


(126) 


But  the  analytic  continuation  of  this  integrand  has  no  singularities  anywhere  in  the  finite  z-plane. 
The  subtraction  procedure  in  equation  (120)  has  accounted  for  all  the  poles  of  N(z)/D(z )  in  the 
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subtracted  term  N(za)/D(z).  Therefore,  numerical  evaluation  of  the  integral  in  equation  (126) 
can  proceed  without  fear  of  encountering  any  singularities  on  or  near  the  path  of  integration. 


At  this  point,  equations  (109),  (125),  (126),  and  (122)  yield 


1 1+7  if  \b\  <  M 

/i“i  if  H<NJ 


(127) 


When  this  information  is  coupled  with  equation  (110),  there  follows  the  major  simplification 


(128) 


which  holds  for  all  relative  sizes  of  \a\  and  \b\.  The  only  cases  where  the  integral  for  /,  in 
equation  (126)  cannot  be  used  is  when  a  =  0  or  b  =  0.  However,  both  of  these  results  are  already 
known  in  closed  form,  namely, 


£>w(0,h)  =  exp 


b2  \  M_-  1 


y  _L 


QmW)  =  1. 


(129) 


The  quantities  needed  for  determination  of  the  integrand  of  equation  (126)  are  given  in 
equations  (111),  (112),  and  (1 19).  The  integral  in  equation  (126)  is  to  be  used  in  conjunction 
with  equation  (128)  when  \b/a\  is  fairly  near  unity.  Otherwise,  the  earlier  form  in  equation  (109) 
is  to  be  used  directly.  The  reason  that  equations  (126)  and  (128)  are  not  used  all  the  time  is  that 
form  (109)  retains  full  significance  even  when  QM(a,b )  and  qM(a,b)  become  extremely  small 
or  large,  whereas  form  (128)  retains  Ml  significance  only  when  QM(a,b )  and  qM{a,b )  are  in 
the  neighborhood  of  0.5.  When  form  (109)  would  tend  to  underflow  or  overflow,  the  logarithm 
of  /  can  be  easily  computed  because  the  leading  exponential  factor  is  the  major  cause  of  that 
problem;  thus,  log (QM)  and  log (qM)  can  be  output  in  these  cases,  in  addition  to  QM  and  qM. 

Atx  =  0  in  equations  (109)  and  (126),  the  exp[-a&(l-cos(x))]  term  is  1.  However,  asx 
increases  from  0  toward  n. ;  this  term  can  grow  very  fast  in  magnitude,  reaching  value  exp(-2ah). 
To  avoid  this  growth,  the  polarity  of  product  a  b  can  be  switched,  if  necessary,  so  that 

real(a  6)  >  0  (130) 

before  equations  (109)  or  (126)  are  numerically  evaluated.  This  switch  is  permissible  according 
to  the  results  in  equation  (104).  The  remaining  quantities  in  the  two  different  integrands  involve 
cos  functions  evaluated  along  the  real  axis.  And  the  dangers  of  a  small  denominator  D(x)  have 
been  avoided  by  employing  equation  (126)  when  |a|  « \b\.  The  only  difficult  case  not  covered  by 
this  switch  is  when  |imag(a  b)\  is  large,  because  the  exp[-oZ>(l  -  cos(x))]  term  will  then 
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oscillate  very  quickly  with  x,  even  as  its  magnitude  decreases.  This  will  necessitate  numerous 
integrand  evaluations  in  order  to  attain  an  accurate  integral  result.  In  the  very  special  case  where 
arguments  a  and  b  are  complex,  but  their  product  is  real,  this  oscillatory  problem  will  not  arise  at 
all;  however,  switch  (130)  is  still  recommended.  The  oscillations  of  the  remaining  cos  functions 
in  the  integrand  of  equations  (109)  and  (126)  are  tolerable  for  moderate  values  of  M. 


Another  case  that  requires  special  treatment  is  when  b  =  a  or  b  =  -a.  Reference  to  equation 
(104)  and  the  closed- form  result  in  reference  4,  equation  (7),  yields 


M- 1 


Qm  ( a>±a )  =  -  +  exp  (-a2 )  -  /„  (a2 )  +  £  Im  (a2 ) 


m= I 


(131) 


Instead  of  subtracting  out  the  constant  N(z0  )  in  the  numerator,  as  done  in  equation  (120),  it 
is  possible  to  subtract  other  factors  instead.  If  so,  a  useful  integral  result,  relative  to  equation 
(109),  is 


cos(M  x)-r  cos ((M  - 1)  x ) 
r  +  — 1  -  cos(x) 


—  rM~x 

for  |r|  >  1 

--1 

for  |r|  =  1 

0 

for  |r|  <  1 

(132) 


A  MATLAB  routine  was  written  to  incorporate  the  above  developments;  it  is  called 
according  to 

[Q,  q,  Q  log,  q  log]  =  Qq  log (M,  a,b) .  (133) 

This  routine  outputs  the  four  quantities  QM,  qM,\og(QM),  and  log  (qM).  At  least  one  of  these 
quantities  always  retains  full  significance.  However,  cases  arise  where  the  two  quantities  that 
have  to  be  subtracted  from  each  other  are  much  larger  than  the  difference,  thereby  losing  all 
significance  in  the  end  result. 
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STEEPEST  DESCENT  EVALUATION  OF  INTEGRALS 


The  basic  problem  of  interest  here  is  to  evaluate  a  contour  integral 

J=j  dz^iz),  (134) 

c 

where  C  is  a  Bromwich  contour  extending  from  -  zoo  to  +  zoo.  It  is  presumed  that  contour  C  can 
be  moved  so  as  to  pass  through  a  (simple)  SP  of  the  integrand  'P(z).  Define 

y/{z)  =  log[vP(z)];  then  'P(z)  =  exp[$/(z)]  =  exp[^r  (z)]  exp[z  y/i (z)],  (135) 

where  subscripts  r  and  i  denote  the  real  and  imaginary  parts,  respectively.  The  locations  {zs }  of 
the  SPs  of  the  integrand  'P(z)  are  located  where 

vXzs)  ~  0-  036) 

Let  Cs  denote  a  steepest  descent  (SD)  contour  passing  through  an  SP  at  zs .  Then,  a  very 
important  property  of  analytic  functions  on  a  path  of  SD  is  that 

Vi(z)  =  Vi(zs)  f°r  2  on  Cs.  (137) 

Thus,  from  equation  (135),  the  angle  ^,  (z)  of  integrand  vF(z)  on  contour  Cs  remains  equal  to 
its  value  at  SP  location  zs.  Equation  (137)  is  the  key  to  finding  the  SD  contours  out  of  an  SP  at 

Define  difference  function 

<fi(z)  =  y/(z)-y/(zs)  for  all  z.  (138) 

Then, 

(f>i (z)  =  Vi (z)  -  y/i (zs )  =  0  for  z  on  C,.  (139) 

An  alternative  way  of  stating  this  property  is 

(f)(z )  is  real  for  z  on  Cs ,  and  ^(zJ)  =  0.  (140) 

Then,  using  equation  (138),  the  integrand  'P(z)  in  equation  (135)  can  be  written  as 


'F(z)  =  exp(>(zs )]  exp[^(z)]. 


(141) 


Substitution  in  equation  (134)  yields 

/  =  |  dz  ¥ (z)  =  |  dz  'P(z)  =exp[^(z5 )]  J  dz  exp[^(z)].  (142) 

c  cs  cs 

On  SD  contour  Cs ,  real  function  <j>{z)  decreases  monotonically  from  the  value  0  as  z  moves 
away  from  SP  location  zs ,  in  either  direction.  Although  integrand  exp[^(z)]  is  real  on  contour 
Cs ,  the  integral  on  z  is  complex  because  dz  is  complex.  The  integrand  of  equation  (142)  at 
z  =  zs  is 

exp0(z,)]  =  exp[O]  =  l,  (143) 

which  is  its  maximum  value  anywhere  on  contour  Cs . 

The  additional  factor  of  exp[y/(zs)]  in  equations  (141)  and  (142)  is  a  complex  constant  and 
can  be  ignored  during  the  evaluation  of  the  contour  integral  in  equation  (142).  However,  this 
factor  may  cause  underflow  or  overflow  in  result  I.  In  such  cases,  the  logarithm  of  I  can  be 
computed  according  to 


log(/)  =  ^(zJ  +  log 


|  dz  exp[^(z)] 


(144) 


Although  the  log  of  the  integral  will  produce  the  standard  principal  value  logarithm,  the 
remaining  term  may  cause  the  imaginary  part  of  log(i)  to  exceed  the  (-ft,  n\  range.  In  this  case, 
the  imaginary  part  of  log(7)  can  be  subject  to  a  mod  operation,  which  will  return  it  to  the  range 
( -7i. ,  tr],  if  desired.  Of  course,  if  and  when  /  is  finally  computed  by  taking  the  exp  of  log(7),  this 
modification  will  be  irrelevant,  since  exp(/  2ftn)  =  \  for  n  integer. 

The  key  property  for  determining  SD  contour  Cs  out  of  SP  location  zs  is,  from  equation 
(139), 

$(z)  =  0  for  z  on  Cs.  (145) 


To  get  started  on  contour  Cs,  consider  the  neighborhood  of  point  zs.  Near  this  SP  location  zv, 
the  following  approximation  is  very  useful: 


<j){z)  =  y/(z)  -  y/(zs )  =  ^  i//"(zs )  (z  -  zv )2 , 


(146) 
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since  y/\zs)  =  0,  and  presuming  y/"(zs)^  0.  For  complex  point  zx  (near  zs)  to  be  on  contour 
Cs,  one  must  have  (in  keeping  with  equation  (145)), 

^\//”(zs)(zx  ~zs)2  =-fs,  fs  real  and  positive.  (147) 

That  is, 


=zs±K 


2fs 


r'M 


(148) 


These  two  values  correspond  to  the  two  directions  of  the  SD  contours  out  of  the  (simple)  SP  at  zs . 
Then,  from  equations  (146)  and  (147), 

VKzi)  =  -/s>  and  exp[^(z,)]  =  exp(-/s).  (149) 

For  example,  if 

fs  =  0.25,  then  exp[^(z,)]  =  exp(-0.25)  =  0.78,  (150) 

versus  the  SP  value  of  exp[^(zs )]  =  1 .  Equation  (148)  can  be  used  to  approximately  find  the  first 
point(s)  on  Cs  away  from  SP  location  zs,  using  fs  of  the  order  of  0.25. 

The  point(s)  zr  given  by  the  right-hand  side  of  equation  (148)  will  not  be  exactly  on  SD 
contour  Cs .  Therefore,  a  local  search  in  the  neighborhood  of  zr  will  be  necessary  to  locate  the 
point  Zj ,  where  ^,  (z, )  =  0.  One  possible  procedure  is  to  let  the  points  (for  the  +i  alternative  in 
equation  (148)) 


z,  =  z„  +  z 


~~  exp(±/^),  ^-^-  =  18  degrees, 

¥  (zs)  10 


(151) 


which,  hopefully,  lie  on  opposite  sides  of  contour  Cs ,  and  compute 

(j){z_)  and  (f>(z+ ).  (152) 

Then,  since  (f),  (z)  =  0  on  Cs,  solve  for  the  linearly  interpolated  value  of  the  zero  crossing  as 


Q  d(Q  +  6(Z+)  Q 


(153) 


Finally,  take  the  first  point  away  from  SP  location  zs  as 
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Z,  =  Z+l 


V"0S) 


exp(/0, ). 


(154) 


Then,  <j>i  (z1 )  =  0  is  expected,  although  not  exactly.  A  similar  procedure  is  used  for  the  —i 
alternative  in  equation  (148),  which  corresponds  to  the  opposite  direction  of  SD  out 
of  the  SP  at  zc. 


Once  the  movement  away  from  the  SP  location  has  taken  place,  a  different  procedure  is  in 
order.  At  two  close  points  zn_l,zn  on  Cs, 


</>(z„ )  =  ^(z„_,  +  [z„  -  z„_,  ])  £  ^(z„_, )  +  (j)\zn_ )  (z„  -  z„_, ). 


(155) 


Presuming  <j){zn_ j)  is  real,  then  to  keep  (j>(zn)  real,  one  must  take 


1)  0„  ~zn-i)  =  -f>  f  real  and  Positive. 


(156) 


That  is, 


z  =  z  t 

n  n-t 


f(v.) 


¥\zn-x) 


(157) 


Then,  the  new  integrand  value  is 

exp[^(z„ )]  £  exp[^(z„_, )  -  /]  =  exp[^(z„  , )]  exp(-/). 


(158) 


For  example,  if  /  =  0.25,  the  relative  amplitude  of  the  integrand  at  the  two  close  points  on  Cs  is 
down  by  exp(-0.25)  =  0.78.  Equation  (157)  cannot  be  used  with  the  SP  location  zs  as  a  starting 
point  because  y/'(zs )  =  0,  and  equation  (155)  is  not  adequate  there.  An  extension  of  equations 
(155)  through  (157)  to  second  order  is  given  in  equations  (166)  and  (167). 

As  point  z„  moves  away  from  SP  location  zs,  ¥'(z„)  possibly  will  become  small, 
suggesting  large  increments  z„  -  zn_ l  according  to  equation  (157).  To  prevent  this,  reconsider 
equation  (148).  The  magnitude  of  the  first  increment  away  from  the  SP  at  zs  is 


z,  -z.  £ 


lku)| 


=  A„. 


(159) 


Then,  a  possible  alternative  to  equation  (157)  is 


z  =z  -A 

n  n— 1  s  t /  \ 


y(z„-x)\  ’ 


(160) 
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which  retains  the  direction  information  of  equation  (157)  and  does  not  allow  large  increments  to 
develop.  In  general,  one  could  use  the  rule 


z„  =  z„- i-g(z„-i)> 


(161) 


where  |g(z„_1)|  has  the  smaller  magnitude  of  the  two  update  terms  in  equations  (157)  and  (160). 
Since  new  point  zn  will  not  lie  exactly  on  SD  contour  Cs ,  let 


z±  =  z«-i  -  g(z„- 1 )  exp (±i0),  6  ~  —  =  1 8  degrees. 


(162) 


Compute  <t>(z±)  and  interpolate  as  in  equations  (152)  and  (153).  Finally,  take 


z„  =z„-i  -g(z„-i)exp(/<9,). 


(163) 


Again,  ^,  (z„)  =  0  is  expected,  although  not  exactly.  This  procedure  requires  two  (f>{z) 
calculations  at  each  stage  (at  z±)  instead  of  one.  However,  it  continually  corrects  itself  by 
maintaining  cf)i (zn )  =  0  for  each  new  contour  point  zn.  This  recursion  procedure  is  repeated 
until  the  real  part  of  <f>{z)  is  sufficiently  negative  that  the  integrand,  exp[^(z)],  in  equation  (142) 
is  negligible  or  below  a  specified  tolerance. 

Once  the  SD  paths  out  of  an  SP  location  have  been  determined,  that  is,  the  sequence  of 
values  {zn}  on  both  sides  of  the  SP,  a  numerical  integration  procedure  is  in  order.  The  sequence 
of  points  {zn }  will  be  close  to  the  path  of  SD,  but  they  will  not  lie  directly  on  it.  Nevertheless, 
analytic  integration  along  the  piecewise  linear  path  afforded  by  locations  (zn }  would  yield  the 
exact  same  result  for  the  integral  as  for  the  SD  path,  according  to  Cauchy’s  theorem.  However, 
since  an  analytic  integration  is  not  possible  in  general,  some  numerical  integration  rule  must  be 
employed  on  the  set  of  complex  locations  { z„ }  and  the  corresponding  complex  function  values, 
{exp[^(zn)]}.  For  example,  the  parabolic  interpolation  and  integration  procedure  for  three 
adjacent  complex  locations  z,,z2,z3  is  given  by 

Z21  =  Z2  ~  Z1  >  Z32  “  Z3  ~  Z2  >  Z31  =  Z3  "  ZI  > 


rt=^(2 z21-z32),  r2  =  — 

Z21  Z21  Z32 


5  r3  ~  (2  ^32  Z2l)’ 

Z  -JO 


(164) 


sum  =  -\rx  f(zx )  +  r2  f(z2 )  +  r3  /(z3 )], 
o 
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where  / ( z )  is  the  function  to  be  integrated.  When  the  three  points  zx,z2, z3  all  lie  on  a  straight 
line  and  are  equally  spaced,  that  is, 

z21=h,  z32=h,  z3l=2h,  h  complex,  (165) 

then  rx  =2  h,r2  =8  h,r3  =  2h,  and  equation  (164)  reduces  to  Simpson’s  rule  with  complex 
function  samples.  However,  the  sequence  {zn }  yielded  by  recursion  (163)  will  generally  lie  on  a 
curved  series  of  points  near  Cs  in  the  z-plane,  necessitating  the  use  of  the  general  rule  in 
equation  (164). 

Since  the  spacing  between  sample  points  {zn }  is  not  zero,  the  numerical  procedure 
employing  equation  (164)  will  not  give  an  exact  result  for  the  desired  integral.  However,  by 
sampling  halfway  between  the  existing  points  and  evaluating  the  integrand  at  these  intermediate 
points,  re-evaluation  of  the  sum  by  means  of  equation  (164)  (at  twice  as  many  points)  will 
improve  the  accuracy  of  the  approximation.  There  is  no  need  to  duplicate  the  function 
evaluations  that  have  already  taken  place;  only  the  function  values  at  the  halfway  points  need  to 
be  evaluated.  This  halving  procedure  can  be  repeated  until  two  consecutive  approximations  for 
•the  integral  of  interest  differ  by  less  than  some  user-specified  tolerance.  The  amount  of  storage 
is  not  excessive,  although  the  number  of  stored  locations  and  function  values  doubles  with  each 
stage  of  this  procedure. 

According  to  equations  (135)  and  (138),  <j>{z)  is  the  difference  of  two  logarithms.  During 
the  tracking  of  the  SD  contour  according  to  equations  (162)  and  (163),  it  is  possible  that  the 
branch  line  of  the  principal  value  logarithm  log(z)  may  be  crossed.  This  will  yield  an  undesired 
jump  of  ±  2  n  in  the  imaginary  part  of  the  logarithm.  This  jump  would  grossly  affect  the 
interpolation  procedure  in  equations  (152)  and  (153),  because  ^,  (z)  must  be  maintained  near 
zero  on  contour  Cs.  Since  the  z  increments  along  contour  Cs  are  small,  according  to  equation 
(163),  the  simplest  way  to  circumvent  this  discontinuity  problem  is  to  always  force  function 
(j){z)  to  return  an  imaginary  part  in  the  range  (-k,  tt\.  If  function  (f>{z)  contains  other  principal 
value  functions,  such  as  sqrt(z),  this  same  correction  must  be  incorporated.  It  is  not  necessary  to 
modify  function  i//(z),  except  possibly  as  noted  under  equation  (144). 

If  the  expansion  of  <j>(z)  is  used  to  second-order,  instead  of  the  first-order  expansion  in 
equation  (155),  that  is, 


^(z„)M(z„_|)  +  ^'(z„_i)(z„ 


there  follows  the  alternative  recursion 


(166) 
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AO 


(167) 


z  =  z 


rt-1 


2/AV,) 


1- 


f(v,)2 


(As  ^"(zn_,)  — >  0,  this  rule  approaches  that  in  equation  (157).) 

When  equation  (136)  is  solved  for  the  SP  locations,  namely, 

1/00  =  0,  (168) 

there  may  be  many  solutions,  depending  on  the  nonlinearity  of  function  y/'(z).  Some  of  these 
SP  locations  may  be  useless;  for  example,  those  SD  paths  that  lead  into  the  locations  of  essential 
singularities  or  zeros  may  not  be  useful  at  all,  or  they  may  have  to  be  augmented  by  an  additional 
path  leading  over  other  SP  locations  on  the  way  to  the  desired  eventual  destination  in  the  z-plane. 
In  other  cases,  it  may  be  necessary  to  use  more  than  one  SP  passage  to  piece  together  a  complete 
path  going  over  the  desired  extent  in  the  z-plane.  These  alternatives  can  be  discovered  only  by 
plotting  the  SD  paths  for  the  various  SPs  and  choosing  the  desired  one(s).  Even  then,  some  of 
the  contour  integral  results  for  /  may  lead  to  obvious  results  and  not  be  useful.  Some 
experimentation  is  required  to  find  the  right  combination(s)  for  the  particular  problem  at  hand. 
However,  as  noted  in  equation  (144),  underflow  and  overflow  can  always  be  kept  under  control, 
and  very  accurate  numerical  results  can  be  obtained. 


EXAMPLE  OF  NON-CENTRAL  CHI-SQUARED  VARIATE 


The  non-central  chi-squared  variate  has  the  PDF  (see  equation  (9)) 


p(u )  =  a 


r  [z — \M~l 

yJ2au 


exp 


■au  — 


P 


2  \ 


1m-\ (P^2au)  for  w  >  0, 


where  a  and  /?  are  real  and  positive.  The  corresponding  MGF  is 


//(A)  =  jdu  exp(Aw)  p(n)  = 


(  3  VM  ( 


1 - 

V  a) 


exp 


A  p 


2  \ 


ya-A  2  j 


for  A,.  =  real(A)  <  a, 


while  the  EDF  is 


(169) 


(170) 


e(v)  =  J  du  p(u)  =  Qm (p, y]2av)  for  v  >  0. 

V 

The  (analytic  continuation  of  the)  MGF  //(A)  has  an  ES  at  A  =  a. 


(171) 
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Instead  of  trying  to  evaluate  the  EDF  by  means  of  integral  (171),  an  alternative  is  available 
by  using  the  inverse  Laplace  transform  result 

p(u)  =  — f  dX  exp(-w  X)  p(X),  ( 1 72) 

Hji  ' 


where  Bromwich  contour  C,  must  initially  be  in  the  region  of  analyticity  of  the  MGF,  namely,  to 
the  left  of  the  ES  at  X  =  a.  Then,  for  v  >  0,  the  EDF  is  given  by 

e(v)=  f  dup(u)=—^—  fdXp(X)  \du  exp(-w/l)  =  — f  —  p(X)  exp(-v^),  (173) 

v  il71 1  v  i2*c\  * 


where  contour  C2  must  now  pass  between  the  pole  at  X  =  0  and  the  ES  at  X  =  a.  The  additional 
restriction  is  required  for  the  latter  u  integral  in  equation  (173)  to  converge,  that  is,  Xr  >  0. 

Now  substitute  equation  (170)  into  equation  (173),  and  equate  the  result  to  equation  (171), 
thereby  getting  an  explicit  integral  relation  for  the  generalized  QM  function  as 


r  dX  ,  . . 

r  xr 

f  X  J32' 

f—  exp(-v  X) 

1 - exp 

d  A 

C2 

l  a) 

^ a-X  2  j 

(174) 


Now  make  the  substitutions 


a  =  — ,  (5  =  a  ( a  and  b  real  and  positive) 
2v 


(175) 


to  get  the  form 


0*  (*»&)  =  -““  f  ~  exp(-v  X) 

1  )  TT  »  A 


H7T  J  X 


1- 


2vX 
b  j 


exp 


X  a 2 


b2/( 2v)-X  2 


(176) 


Finally,  let  X  -  z  —  and  get 
2v 


0w  (*>*)  =  7^;  jy  exp(-5z)(l-z)“M  exp^ 


Az  ' 

1  ~Zj 


(177) 


where 

A  =  a2/2,  B  =  b2/2, 


(178) 
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and  contour  C3  must  pass  between  the  pole  at  z  =  0  and  the  ES  at  z  =  1  of  the  integrand  of 
equation  (177).  Under  the  identification 


1  (  At  ^ 

^(z)  =  —  exp(-i?z)  (1-z)  M  exp  - 

z  W-ZJ 


the  integral  of  interest  takes  the  form 
Q«(a,b)  =  ~I,  /=f<fe4>( z). 

i2n  “ 


(179) 


(180) 


The  residue  of  T'(z)  atz  =  0  is  obviously  1,  as  seen  from  equation  (179);  that  is,  Res(0)  =  1. 
Since  B  is  real  and  positive,  the  contour  C3  can  be  moved  far  into  the  right-half  z-plane,  where 
the  integrand  vP(z)  tends  to  zero.  The  only  singularity  encountered  in  this  movement  is  the  ES 
at  z  =  1 .  Therefore,  I  =  -i  2  n  Res(l),  where  the  minus  sign  is  due  to  the  negative  (clockwise) 
sense  of  the  encirclement  of  the  ES  atz  =  1.  Thus,  equation  (180)  yields  the  very  useful 
observation  that 


Res(l)  =  -0A/(fl,ft).  (Also,  Res(0)  =  l.)  (181) 

The  situation  can  now  be  summarized  as  follows.  Any  contour  C  that  encircles  just  the  pole 
at  z  =  0  (in  a  positive  sense)  will  yield  the  value  1  for  the  integral 

J  =  —  f  dz  'F(z).  (182) 

i2n  * 

Also,  any  contour  C  that  encircles  just  the  ES  at  z  =  1  will  yield  the  value  -  QM  (a,  b )  for  J.  Any 
contour  C  that  encircles  both  singularities  will  yield  the  value  1  -  QM  (a,  b )  for  J.  These  observa¬ 
tions  will  allow  easier  interpretation  of  the  numerical  results  for  the  SPs  and  associated  SD  paths 
to  be  presented  in  the  examples  below. 

According  to  equations  (135)  and  (179),  the  function 

Viz)  =  logPF(z)]  =  -Bz  +  -  log(z)  -  M  log(l  -  z),  (1 83) 

1  -  z 


where  A  and  B  are  given  by  equation  (178).  This  function  in  equation  (183)  uses  two  principal 
value  logarithms;  therefore,  the  correction  procedure  on  <j>{z)  =  y/(z)  -  i//(zs ),  outlined  two 
paragraphs  below  equation  (165),  must  be  employed.  Then, 
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(184) 


B  zi  +(M  +  \-2B)z2  +(B-A-2-M)z  +  \ 
z(l-z)2 


Saddlepoint  locations  {zs}  of  integrand  'P(z)  are  determined  by  the  three  solutions  of  the  cubic 
equation 

Bz]  +  (M  +  l-2B)z]  +(B-A-2-M)zs  +1  =  0.  (185) 

Although  this  cubic  equation  has  an  analytic  solution,  it  is  not  a  simple  one,  even  for  low-order 
integers  M.  Rather,  given  numerical  values  of  A,  B,  and  M,  the  three  SP  locations  {zs }  can  be 
determined  numerically.  It  is  then  necessary  to  determine  the  SD  contours  leading  out  of  each  of 
these  SP  locations,  and  to  determine  which  of  these  contours  are  useful,  according  to  the 
discussion  accompanying  equation  (182). 

For  A  and  B  real  and  positive,  it  is  seen  from  equation  (179)  that,  on  the  real  axis,  vP(z)  is 
negative  to  the  left  of  the  origin  and  has  a  single  maximum  in  that  region;  that  is,  its  magnitude 
has  a  single  minimum.  In  the  region  between  z  =  0  and  z  =  1,  'P(z)  is  positive  and  has  a  single 
minimum.  To  the  right  of  z  =  1,  the  magnitude  of  'P(z)  has  a  single  maximum ;  observe  that  the 
value  of  the  ES  term  rapidly  approaches  zero  as  z  -» 1  + .  These  three  locations  correspond  to 
the  SP  locations  for  this  case.  However,  the  SD  directions  for  the  third  SP  are  along  the  real  axis 
of  the  z-plane  and  one  path  heads  directly  into  the  ES  at  z  =  1,  making  this  SP  useless. 

For  the  SP  located  between  0  and  1,  on  the  other  hand,  the  SD  paths  start  out  perpendicular 
to  the  real  axis,  and  then  both  bend  to  the  right,  due  to  the  exp  {-Bz)  term.  This  is  contour  C3 
in  equation  (180).  The  only  singularity  encircled  (in  the  negative  sense)  by  this  contour  is  that  at 
z  =  1,  yielding,  upon  use  of  equation  (181), 

\dz  T(z)  =  -Res(l)  =Qm (a,b),  (186) 

i2n  “ 

which  is  consistent  with  equations  (180)  and  (182). 

For  the  SP  location  to  the  left  of  the  origin,  the  SD  paths  also  start  out  perpendicular  to  the 
real  axis,  and  then  both  bend  to  the  right,  again  due  to  the  exp(-2?z)  term.  Call  this  contour  C4. 
Therefore,  both  singularities  of  'P(z)  are  now  encircled  (in  the  negative  sense),  leading  to  a 
value  for  the  integral  of 

TT  f  dz  ¥(z)  =  -Res(0) - Res(l)  =  -1  -[~Qu (a,b)]  =  QM(a,b )  - 1.  (187) 

i2n  • 


y/'(z)  =  -B  + 


(1-z) 


2  z  1-z 


2A  1  M 

¥  (z)  ~ - r  +  — r  + - 7 

(1-z)3  z2  (1-z)2 


37 


This  is  consistent  with  equation  (182).  Thus,  contour  C4  through  the  leftmost  SP  location  yields 
not  Qm (a, 6),  but  rather  QM (a, b)-l.  This  can  be  a  very  useful  observation  when  QM (a, b)  is 
very  close  to  1 .  The  complementary  function 


qM(a>b)  =  l-QAf(a’b)  O88) 

is  the  cumulative  distribution  and  can  be  evaluated  very  accurately  far  out  on  its  tail  values  by 
means  of  contour  C4,  whereas  its  evaluation  by  the  difference  in  equation  (188),  coupled  with 
contour  integral  (186)  on  contour  C3,  could  have  no  significance  at  all.  This  example  illustrates 
that  some  investigation  is  in  order  for  the  various  SPs  and  their  associated  paths  of  SD,  in  order 
to  find  out  which  is  the  most  useful  path  for  the  problem  at  hand. 


NUMERICAL  EXAMPLES 

A  series  of  numerical  examples  of  equation  (179)  is  now  presented  that  will  further 
demonstrate  the  utility  and  futility  of  some  of  the  various  SPs.  The  first  case  is 

M  =  3,  a  =  1.1,  6  =  2.1;  zs  =-1.2866,  0.30086,  1.1716.  (189) 

This  case  bears  out  all  the  observations  made  above.  The  SP  locations  {zs }  will  be  labeled  1 ,2,3 
in  the  order  presented  here.  Also,  Qu(a,b )  will  be  abbreviated  by  Q  in  these  examples. 

The  second  case  is 

M  =  3,  a  =  1.1,  6  =  21;  zs  =-0.00461,  0.94054,  1.0459.  (190) 

Since  6  is  much  larger  then  a,  there  follows  Q  =  3.001  e  —  85  —  i  1 .352e  — 101,  which  was 
obtained  by  integration  along  the  path  passing  through  the  second  SP  (at  zs  =  0.94054 ).  This 
path  encircled  the  ES  in  the  negative  sense,  thereby  yielding  equation  (186)  again.  The  path 
through  the  first  SP  encircled  both  singularities  in  the  negative  sense,  and  should  yield 
-Res(0)-Res(l)  =-l  +  Q,  which  is  virtually -1.  The  numerical  integration  confirmed  this 
result;  however,  this  is  not  a  useful  path  or  result  due  to  the  loss  of  significance.  The  third  SP 
was  useless,  with  an  SD  path  leading  into  the  ES. 

The  third  case  is 

M  =  3,  a  =  21,  6  =  1.1;  ^=-21.725,0.00445,17.109.  (191) 

Since  a  is  much  larger  than  6,  there  follows  -q  =  -1  +  Q  =  — 1 .047e  —  9 1  —  i  3. 1 98e  —  1 07,  which 
was  obtained  by  integration  along  the  path  passing  through  the  first  SP.  This  path  encircled  both 
singularities  in  the  negative  sense,  yielding  -  Res(0)  -  Res(l)  =  -1  +  Q  =  -q,  as  expected.  The 
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path  through  the  second  SP  encircled  only  the  pole  at  the  origin  in  a  positive  sense,  yielding 
integral  value  Res(0)=l.  This  is  correct  but  not  useful.  The  third  SP  was  useless. 


The  expression  for  QM(a,b)  in  equations  (177)  and  (178)  is  analytic  in  a  and  b  for  all  a  and 
b.  Therefore,  by  use  of  analytic  continuation,  the  integral  expression  can  be  used  for  complex 
values  of  a  and  b.  Accordingly,  the  fourth  case  is 

M  =  3,  a  =  1.1,  b  =  2.1  +  /  2; 

(192) 

z,  =0.1425  +  / 0.1240,  0.6374  +  i 0.8624, 1.1738-/0.0362.  7 

Since  b  is  complex,  the  SP  locations  can  now  lie  in  the  complex  z-plane,  not  just  on  the  real  axis. 
The  path  through  the  first  SP  starts  and  ends  in  the  fourth  quadrant  and  encircles  only  the  pole  at 
the  origin  in  a  positive  sense;  the  value  of  the  integral  is  1,  which  is  not  a  useful  result.  The  path 
through  the  second  SP  encircles  both  singularities  in  the  negative  sense,  yielding  integral  value 
Q  - 1  =  -q  =  2. 1 123  -  i  4.6039.  The  path  through  the  third  SP  heads  into  the  ES  and  is  not  useful 

The  fifth  case  is 


M-  3,  a  =  l.l  +  i  3,  b-2.\/a-, 

zs  =0.0495-/0.2181,  0.2505  +  i  0.9875,  15.831  -i  12.742. 


(193) 


Here,  both  a  and  b  are  complex,  although  their  product  is  real.  The  path  through  the  third  SP 
encircles  both  singularities  in  a  positive  sense,  yielding  value  Res(0)  +  Res(l)  =  1  -  Q  =  q,  which 
is  a  useful  result.  The  path  through  the  first  SP  starts  in  the  second  quadrant  and  ends  up  heading 
into  the  ES  at  z  =  1;  by  itself,  this  is  not  useful.  However,  the  path  through  the  second  SP  comes 
out  of  the  ES  and  heads  off  into  the  second  quadrant.  The  combination  of  these  two  paths 
encircles  the  pole  in  the  positive  sense,  but  does  not  encircle  the  ES,  though  these  two  SD  paths 
both  enter  into  the  ES  at  the  same  angle.  Therefore,  the  sum  of  the  two  integrations  obtained  by 
using  the  combination  path  is  equal  to  Res(0)=l .  This  is  an  interesting  combination  of  SD  paths, 
but  it  does  not  lead  to  a  useful  result. 

The  sixth  case  is 


M  =  3,  a  =  10,  b  =  H0;  zs  =  0.009615,  1.0352 ±i  1.0042.  (194) 

The  path  through  the  first  SP  encircles  the  origin,  yielding  value  Res(0)=l.  As  above,  the 
combination  of  the  two  paths  through  the  two  remaining  conjugate  SPs  encircles  both 
singularities  in  the  positive  sense,  yielding  Res(0)  +  Res(l)  =  1  -  Q  =  q  =  -0.026941 .  The  path 
through  the  upper  SP  stays  in  the  upper-half  z-plane,  while  the  path  through  the  lower  SP  stays  in 
the  lower-half  z-plane. 
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The  seventh  case  is 


M-  3,  a  =  10  +  /,  b  =  0.5  +  / 10; 

z,  =0.00965-/0.00048,  0.8848  +  / 1.001,1.1849-/ 0.9926. 


(195) 


Now,  both  parameters  a  and  b  are  complex,  and  the  SP  locations  have  no  symmetry  properties. 
The  path  through  the  first  SP  encircles  just  the  origin,  yielding  integral  value  1  as  usual.  The 
path  through  the  third  SP  starts  in  the  fourth  quadrant  and  works  its  way  toward  the  positive  real 
axis;  it  ends  up  in  the  ES  location,  but  the  function  values  get  so  small  that  computation  ceases  at 
the  specified  tolerance  before  that  could  happen.  The  path  through  the  second  SP  starts  near  the 
positive  real  axis  and  ends  up  in  the  second  quadrant;  it  starts  at  the  ES  location,  but  the  function 
values  are  again  too  small.  Each  of  these  last  two  contours,  by  themselves,  do  not  yield  a  useful 
result.  However,  the  combination  of  the  two  paths  encircles  both  singularities  in  the  positive 
sense,  leading  to  integral  value  Res(0)  +  Res(l)  =  1  -  Q  =  q  =  1.7196-  /  5.2731.  Thus,  here  is  an 
example  where  a  combination  of  two  SD  paths  must  be  used  to  attain  a  useful  result.  Taken 
separately,  none  of  the  three  SD  contours  through  the  three  SPs  yield  useful  integral  results 
relative  to  QM(a,b )  or  qM(a,b).  In  fact,  it  is  conceivable  that  in  some  cases,  it  may  be 
necessary  to  use  a  combination  of  all  three  SPs  and  their  respective  descent  paths  to  achieve  a 
meaningful  integral  result. 

The  eighth  case  is 

M  =  1000,  a  =  1 . 1  +  /  3,  b  =  2A-i  2, 

(196 

zs  =-10.602-/237.77,  0.00100-/ 7.5455e- 6,  0.9961 1  +  / 0.0033. 


The  SD  paths  through  the  first  SP  encircle  both  singularities  in  a  negative  sense.  However, 
y/{zs )  =  -4473.2  +  / 1.0377  has  such  a  large  negative  real  part  that  the  factor  exp[^/(zs)J  in 
equation  (142)  underflows.  In  this  case,  it  is  necessary  to  resort  to  the  procedure  mentioned  in 
equation  (144);  the  end  result  is  that  log(7)  =  -4472.2  +  /  2.5598,  which  retains  significance  for 
further  computations.  Use  of  the  second  SP  leads  to  encirclement  of  just  the  pole  at  the  origin 
and  a  known  value  of  1  for  the  integral.  The  path  associated  with  the  third  SP  leads  into  the  ES 
location  and  is  not  useful.  This  example  illustrates  that  large  values  of  Mean  be  handled  by  this 
combined  SP  and  SD  procedure,  although  it  may  be  necessary  to  output  a  logarithm  of  the 
integral  instead  of  the  integral  value  itself. 

The  basic  difficulty  with  this  approach  is  the  presence  of  the  ES.  The  ES  acts  like  a 
directive  “black  hole,”  attracting  any  nearby  SD  path  directly  into  itself.  The  SD  paths  through 
the  three  SPs  are  adversely  affected  by  this  effect,  making  interpretation  of  these  results  very 
difficult  and  virtually  impossible  to  automate.  Visual  inspection  of  all  the  SD  paths  is  required 
in  order  to  make  the  correct  decision  about  what  path(s)  to  use.  Integrand  (179)  has  no  zeros 
anywhere  in  the  finite  z-plane;  therefore,  the  SD  paths  cannot  terminate  in  zeros,  fortunately. 
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AN  ALTERNATIVE  STEEPEST-DESCENT  APPROACH  FOR  Q^b) 


The  presence  of  both  a  pole  and  an  ES  in  the  integrand  'P(z)  in  equation  (179)  causes 
difficulties  in  interpretation  of  the  SPs  and  their  associated  paths  of  SD.  It  is  desirable  to  have  an 
alternative  integrable  form  that  does  not  have  the  ES.  Such  a  form  can  be  obtained  as  follows. 

In  the  top  line  of  equation  (108),  let  y  =  - x  in  the  latter  half  of  that  equation.  Then,  for 
a  ^  0,  b  ^  0,  \a\  <  |Z>|,  the  resultant  two  halves  of  the  equation  can  be  combined  into  the  form 


QM(a’b)=  exP 


f  a2+b2')fb'\M~l 

2  ){a, 


=  exp 


-Ua-b) 

.  2  j 


b 


i-  fAexpI^cosWl^Pl^-Ja 
17  I  -  exp(~/  a  ) 

exp[-  a  b  [1  -  cos(x)]  +  i  x  ( M  - 1)] 


(197) 


—  f  dx 
In  1 


a 


1  —  —  exp(-z  x) 


Define  the  analytic  continuation  of  the  numerator  of  the  integrand  as 
N(z)  =  exp  [-a  b  [1  -cos(z)]  +  /z  (M  - 1)],  (198) 

which  is  an  entire  function  with  no  zeros  in  the  finite  z-plane.  Also,  define  the  denominator 


D(z )  =  1  -  —  exp(-i  z).  (199) 

b 

Observe  that  N( 0)  =  1,  D( 0)  =  1  -a/b,  and 

N(z  +  2  nn)  =  N(z),  D(z  +  2  nn)  =  D(z)  for  n  =  0,  ±  1,  ±  2, . . . ;  (200) 

that  is,  both  functions  have  period  In  in  the  direction  of  the  x-axis.  Then,  the  integrand  of 
equation  (197)  is 

T(z)  =  ^.  (201) 

D(z) 

There  follows 


( 


a 


i//(z)  =  Log  'P(z)  =  -a  b  [1  -  cos(z)]  +  iz  ( M  - 1)  -  log^l  exp(-iz) 


+  i2nn. 


(202) 
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where  the  two  functions,  Log(  )  and  log( ),  are  the  general  logarithm  and  principal  value 
logarithm,  respectively,  while  n  is  an  arbitrary  integer.  Also, 


y'(z)  =  -a  b  sin(z)  +  i(M- 1)  -  /  -aLbZ*?t'Ji  . 

\-a/b  exp(-iz) 

=  -a  b  sin(z)  +  i  ( M  - 1)  -  i  - — — — 

D(z) 

=  -a  b  sin(z)  +  i  M - - — , 

D(Zy 

y/\z)  =  -a  b  cos (z)  + 


(203) 


(204) 


The  SP  locations  { zs }  of  the  integrand  'P(z)  in  equation  (201)  are  at  the  solutions  of  the 
equation  i//'(zs)  =  0.  Let  e  =  exp (~izs).  Then,  from  equation  (203), 


1/e  -  e  i 

- a  b — - +  iM - 

i  2  \-(a/b)e 


=  0, 


(205) 


which  can  be  simplified  to 

a2  be 3  -a(2M  +  b2)e2  +b  (2M-2-a2)  e  +  a  b2  =0. 


(206) 


Although  this  cubic  in  e  can  be  solved  analytically,  the  solution  is  very  complicated,  even  for 
small  values  of  integer  M.  Instead,  solve  this  cubic  numerically  for  the  three  values  of  e.  Then, 

zs  =  i  Log(e)  =  /  log(e)  +  2nn,  n  -  0,±  1,±  2,...,  (207) 


gives  the  locations  of  the  three  SPs  in  each  vertical  strip  of  width  2  n  in  the  z-plane. 
The  poles  of  integrand  'F(z)  are  located  where  D(z  )  =  0;  from  equation  (199), 


j  exp (-/ zp)  =  1,  exp(i zp )  =  ^ ,  izp=  Log^  =  log| 


+  i2ftn, 


(a\ 

z  =2^r«  +  angle  —  -/log 

\b) 


a 


(208) 


This  pole  location  is  independent  of  M.  There  is  only  one  pole  in  each  2  n  vertical  strip.  And 
since  |a|  <  |h|,  this  pole  location  lies  above  the  real  axis  of  the  z-plane.  The  residue  of  ^(z)  at 
this  pole  location  zp  is 
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N(zp) 

D'(zd) 


exp 


,  1 

(a  b > 

-  a  b  +  a  b  — 

- J - 

2 

\b  aj 

UJ 

.  M-i 


.  a  b 

1 - 

b  a 


(l,  .,0 

tUJ 


(209) 


Therefore,  from  equation  (197), 

V|"_1 

\a; 


exp 


1  ,  „2}(b 


--(a~b) 

v  2  j 


—  f  dz 
2  n  l 


Hz) 

D(z) 


=  1 


(210) 


if  contour  C  encloses  just  the  pole  at  z  ,  in  a  positive  sense. 


LOCATIONS  OF  THE  HILLS  AND  VALLEYS  OF  INTEGRAND  T(z) 

As  y  -»  ±oo,  the  dominant  quantity  in  N(z)  in  equation  (198)  is 
g(z)  s  exp[a  b  cos(z)], 

as  will  be  seen  below.  This  function  has  period  2  n.  Let 
ab  =  n  exp(z'  /3), 

where  n  is  positive  real  and  P  is  real.  Then, 


g(z) =  exp 
=  expj 


exp(z  P)  |exp(/  z)  +  exp(-z  z)} 


exp(z  p  +  ix-y) 


exp 


exp(z p~ix  +  y) 


As  y  — >  +oo,  then 


g(z) ~ exp 


y  exp  [i(p-x)]  exp(y) 


(211) 


(212) 


(213) 


(214) 


which  dominates  exp[z  z  (M  - 1)]  =  exp[(z  x  -  y)(M  - 1)].  Therefore,  the  centers  of  the  hills  of 
g(z)  at  +  zoo  (the  location  of  fastest  growth  of  the  magnitude)  are  at  x  coordinates 

xH+  =  p  =  Angle(a  b)  =  angle(a  b)  +  2n;n,  (215) 
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where  the  two  functions,  Angle( )  and  angle( ),  are  the  general  angle  and  principal-value  angle, 
respectively.  Also,  the  centers  of  the  valleys  of  g(z)  at  +  /oo  (the  location  of  fastest  decay  of  the 
magnitude)  are  at 

xy+  =  (5  +  n  =  angle(a  b)  +  n  +  2  n  n.  (216) 

On  the  other  hand,  as  y  — >  -oo,  then 


g(z)  ~  exp 


£ 

2 


exp[/(/?  +  x)]  exp(-y) 


Now,  the  centers  of  the  hills  at  -  i  oo  are  at 

xH_  =  -  ft  =  -Angle(a  b)  =  -angl e(a  b)  +  2nn , 
while  the  centers  of  the  valleys  at  -  /  oo  are  at 
xv_  =  n  -  P  =  -angle(«  6)  +  #  +  2  #  n. 


(217) 


(218) 


(219) 


MOVEMENT  OF  CONTOUR  FOR  <2^,*) 

The  original  contour  of  integration  for  QM(a,b)  in  equation  (197)  is  a  finite  contour  along 
the  real  axis  of  the  z-plane,  from  -  n  to  n.  However,  since  the  integrand  of  equation  (197)  has 
period  2  n,  the  range  of  integration  can  be  taken  as  xc  to  xc  +  2  n,  where 

xc  =  angle(a  b)  -  n\  (220) 


see  equation  (216).  Then,  two  vertical  lines  can  be  drawn  from  xc  and  xc  +2  k  toward  +  /  oo, 
right  into  the  centers  of  two  adjacent  valleys.  The  combination  of  these  two  vertical  lines  and  the 
portion  of  the  real  axis  between  them  is  a  U-shaped  contour,  to  be  denoted  by  Cu .  Integration 
down  the  left  vertical  line  exactly  cancels  integration  up  the  right  vertical  line,  due  to  the 
periodicity  of  the  integrand.  Thus,  there  follows,  from  equations  (197)  through  (201), 


(  1  , 

QM(a,b)  =  Qxp  --(< a-b ) 

V  z 


b 

yaj 


—  f  dz 
2n  ) 


N{z) 

D(z) 


(221) 


Now,  infinite  contour  Cv  can  be  moved  upward  so  as  to  pass  through  a  convenient  SP(s); 
call  this  new  contour  Cs.  If  the  pole  of  the  integrand  at  zp  (which  is  in  the  upper-half  z-plane 
for  the  present  case  of  |aj  <  |6|)  is  crossed  during  this  contour  movement,  then  use  of  equation 
(210)  reveals  that 
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This  result  is  useful  for  evaluating  complementary  function  qM  (a,  b)  =  1  -  QM  ( a ,  b ).  However,  if 
the  pole  at  z  is  not  crossed  when  moving  contour  Cv  into  Cs ,  then  the  additive  1  on  the  right- 
hand  side  of  equation  (222)  must  be  removed,  and  a  direct  integral  relation  for  Qu  ( a,b )  results. 

In  some  cases,  contour  Cs  may  not  look  like  a  U-shaped  curve.  It  may  be  necessary  to  use 
two  SPs  and  piece  together  two  steepest-descent  contours  that  enter  a  valley  at  -  too  with  an 
x-coordinate  midway  between  the  locations  of  the  original  two  adjacent  valleys  of  Cv  at  +  ice. 

Of  course,  contour  Cs  would  still  start  in  the  left  valley  at  xc  +i  oo  and  end  up  in  the  right  valley 
at  xc  +  2  n  +  i  oo.  Whether  equation  (222),  with  or  without  the  additive  1 ,  should  be  used 
depends  on  the  exact  positions  of  the  three  SPs  and  the  pole  of  integrand  'P(z). 

The  development  above  presumed  that  \a\  <  \b\,  meaning  that  pole  location  zp  lies  in  the 
upper-half  z-plane.  Keeping  a  and  b  fixed  for  the  moment,  contour  Cu  in  equation  (221)  can  be 
moved  downward  in  the  z-plane,  according  to  Cauchy’s  theorem,  without  changing  the  value  of 
the  integral  (221)  for  QM(a,b );  call  this  new  contour  Cd.  Then,  since  QM(a,b )  is  analytic  in 
both  a  and  b,  the  values  of  these  two  parameters  can  now  be  varied  as  desired,  provided  that  the 
pole  location  does  not  cross  Cd\  this  can  always  be  accomplished  by  sufficient  prior  movement 
of  Cd .  Thus,  if  contour  Cd  remains  below  the  pole  location,  the  expression 

n  M  f  1,  M2 1  f  .  N(z) 

Qu(a’h)  =  QxP  -~(a-b)  -  —  I  dz——  (223) 

V  2  )\a)  2k  c*  D{z) 

holds  true  for  all  a  and  b,  not  just  \a\  <  \b\.  Contour  Cd  can  now  be  moved  so  as  to  pass  through  an 
SP(s),  resulting  in  a  modified  contour  Cs .  Again,  relation  (222)  remains  valid  if  the  pole  location 
is  crossed  while  Cd  is  moved  into  Cv ,  whereas  the  additive  1  on  the  right-hand  side  is  dropped  if 
the  pole  location  is  not  crossed  during  movement.  Thus,  relation  (222)  can  be  used  for  all  values 
of  a  and  b,  provided  that  this  “crossing  rule”  is  adhered  to.  In  no  case  does  contour  Cs  in  equation 
(222)  pass  through  the  pole  location  zp ,  because  SD  paths  always  automatically  avoid  poles. 


SELECTION  OF  AN  APPROPRIATE  SADDLEPOINT 

Equations  (206)  and  (207)  give  the  locations  of  the  three  SPs  of  integrand  T'(z)  of  equation 
(201),  while  equation  (208)  gives  the  pole  locations.  When  the  SD  paths  through  the  three  SPs 
are  plotted,  a  wide  variety  of  contours  is  obtained,  making  it  very  difficult  to  decide  what 
combination  to  select,  without  making  a  detailed  observation  of  the  particular  paths.  This  makes 
automation  of  the  Qu(a,b )  evaluation  very  difficult. 
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Typically,  one  SP  is  associated  with  the  pole  and  always  lies  in  close  proximity  to  the  pole 
location.  The  SD  path  through  this  particular  SP  surrounds  the  pole  and  frequently  comes  and 
goes  into  the  same  valley  at  infinity.  Thus,  the  integral  along  this  SD  path  always  gives  the 
values  ±  1  and  is,  therefore,  useless.  The  SD  paths  through  the  two  remaining  SPs  sometimes 
have  to  be  combined,  to  yield  a  complete  path  leading  between  two  adjacent  valleys.  Yet,  on 
other  occasions,  one  SP  will  connect  two  valleys  at  z  =  ±i  oo,  while  the  remaining  SP  will 
accomplish  the  desired  task  of  connecting  two  adjacent  valleys  at  z  =  +z'co.  These  complications 
require  that  an  alternative  tack  be  taken. 

From  this  point  on,  the  complex  parameters  a  and  b  of  QM  (a,  b)  will  be  restricted  so  that 
their  product  is  real.  Then,  equation  (216)  reveals  that  there  are  valleys  at  +  z'oo  located  at  x- 
coordinates  ±  tc.  Instead  of  trying  to  connect  these  two  valleys  at  z  =  +i  oo  by  means  of  an  SD 
path  through  an  SP(s),  a  simpler  approach  will  be  adopted.  Namely,  a  horizontal  path  of 
integration  (not  the  real  axis)  through  an  “appropriate”  SP  over  the  real  interval  ( ~n,n )  will  be 
employed;  the  vertical  legs  of  the  paths  reaching  to  z  =  +i  co  can  be  dropped,  because  the  integrals 
along  these  paths  cancel  each  other.  The  appropriate  SP  is  that  one  with  the  minimum  magnitude 
of  the  slope  of  the  SD  path  as  it  passes  through  the  SP.  This  automatically  eliminates  those  SPs 
that  connect  two  different  valleys  at  z  =  ±/oo,  because  the  slope  for  those  SPs  is  nearly  vertical. 
(Other  alternatives  were  tried,  but  none  proved  to  be  unambiguous  for  all  the  cases  tried.)  Of 
course,  this  horizontal  path  of  integration  is  no  longer  the  path  of  SD  through  that  SP,  but  it  is  still 
a  path  of  descent.  The  peak  magnitude  of  the  integrand  on  that  horizontal  path  occurs  right  at  the 
SP  location,  which  is  a  very  desirable  condition.  This  horizontal  path  through  the  SP  does  not 
eliminate  oscillations  in  the  real  and  imaginary  parts  of  the  integrand,  but  it  does  tend  to  minimize 
oscillations.  This  switch  in  paths  from  SD  to  a  horizontal  path  also  eliminates  the  considerable 
complex  numerical  effort  associated  with  determining  the  actual  SD  path  through  a  particular  SP. 

Since  this  horizontal  path  of  integration  might  intersect  or  come  close  to  the  pole  location, 
indent  the  path  with  a  portion  of  a  small  circle  centered  at  the  pole  location,  where  and  if 
necessary,  thereby  avoiding  a  close  approach  of  the  integration  path  to  the  pole.  Then,  use  a 
numerical  integration  procedure  on  the  horizontal  path  (and  another  numerical  integration 
procedure  on  the  indented  circular  portion  if  present)  such  as  MATLAB’s  quadf  routine.  This 
latter  routine  is  an  adaptive  Gauss-Lobatto  integration  procedure  with  high  accuracy  and 
efficiency  for  smooth  integrands.  By  breaking  the  integration  path  into  two  pieces  (when 
necessary),  the  integrands  in  the  two  quadt?  applications  can  be  kept  as  smooth  functions. 

This  procedure  uses  all  the  information  about  the  locations  of  the  three  SPs  and  the  pole,  but 
avoids  having  to  evaluate  the  detailed  SD  paths  from  any  of  the  SPs.  There  is  no  analytical 
solution  for  the  SD  paths  through  any  of  the  SPs.  Examples  of  the  integrand  on  the  horizontal 
path  reveal  that  moderate  oscillations  are  encountered  as  the  magnitude  of  the  integrand  decreases 
from  its  peak  value  at  the  SP;  however,  severe  oscillations  that  would  cause  cancellations  and  loss 
of  significance  are  not  encountered.  Also,  underflow  and  overflow  can  be  controlled  by  resorting 
to  the  factorization  employed  in  equations  (142)  through  (144). 
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DETERMINING  IF  AN  SD  PATH  PASSES  UNDER  OR  OVER  A  POLE  LOCATION 


It  has  been  seen  above  that  a  key  issue  in  evaluating  the  contour  integral  is  whether  a 
particular  SD  path  passes  under  or  over  a  pole  location  in  the  complex  z-plane.  This 
determination  can  be  accomplished  without  drawing  all  the  various  contours  and  can  be 
automated  as  follows.  Let  the  known  locations  of  all  the  SPs  be  denoted  by  {zs},  s  =  1:  S. 

Also,  let  Cs  denote  the  5-th  SD  contour,  which  passes  through  zs.  Then,  from  equation  (137), 

y/j (z)  =  y/. (zs )  forz  on  Cs,  5  =  1:5.  (224) 

The  pole  location,  zp-xp  +  iyp ,  must  be  known;  for  example,  see  equation  (208).  In  order  to 
determine  if  the  SD  contour  Cs  passes  under  or  over  this  location,  set 

V'i(xp+iyPs)  =  V'i(zs)  (225) 

and  solve  for  yps .  This  solution  must  be  conducted  for  each  5  value  of  interest.  If  yps  >  0,  then 
SD  contour  Cs  passes  above  the  pole  location.  If  yps  <  0,  then  Cs  passes  below  the  pole 
location.  Solution  y  can  never  equal  zero  because  SD  paths  automatically  avoid  any  poles. 

If  there  are  multiple  real  solutions  for  yps  in  equation  (225),  it  will  be  necessary  to  trace  the 
detailed  SD  paths  out  of  the  SP  of  interest  to  determine  exactly  where  the  SD  paths  go.  The  large 
variation  of  path  possibilities  makes  this  computation  unavoidable. 

In  terms  of  the  auxiliary  function  </>{z)  defined  in  equation  (138),  it  is  necessary  to  introduce 


a  more  explicit  notation,  namely, 

<t>siz)  =  y/(z)-y/(zs)  for  5  =  1 :  5.  (226) 

Then,  the  imaginary  parts  become 

4  (z)  =  ¥i  (z)  -  ¥i  (zs )  for  5  =  1 :  5,  (227) 

and,  in  particular, 

</>s,  (xp  +  i  yPs )  =  Vi  (xp  + *  y,, )  -  Vi  (zs )  for  5  =  1 : 5.  (228) 

Therefore,  solution  of  equation  (225)  is  equivalent  to  solving 

UxP+iyPs)  =  0  (229) 

for  yps  for  those  values  of  s  of  interest. 
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SUMMARY 


A  procedure  for  evaluation  of  the  generalized  QM(a,b )  function,  with  complex  arguments 
a,b ,  has  been  derived,  which  keeps  underflow  and  overflow  in  check.  A  MATLAB  routine  is 
listed  in  the  appendix.  It  is  limited  to  the  case  where  the  product  a  b  of  its  arguments  is  real. 
However,  the  procedures  detailed  in  this  report  could  be  used  to  extend  the  routine  to  general 
complex  arguments. 

The  complementary  function  qM  ( a ,  b)  =  1  -  QM  (a,  b)  often  occurs  in  calculations  for  a 
practical  problem.  Accordingly,  it  is  necessary  to  also  calculate  and  supply  this  quantity.  In 
fact,  in  general,  the  four  quantities 

Qm>  Vm*  log (QmX  lo g(qM)  (230) 

should  all  be  output  from  the  general  routine.  This  guarantees  that  at  least  one  of  the  four 
complex  quantities  retains  significance  for  all  argument  values.  This  situation  is  similar  to  the 
case  in  MATLAB,  where  routines  for  both  the  error  function  erf( )  and  the  complementary  error 
function  erfc( )  are  furnished.  All  of  these  features  are  included  in  the  final  routine  for  the 
Qm  {a,  b)  function.  In  particular,  the  call  is 

[Q’<1>Q l°g> <7 log,  A n]  =  Qq(M,a,b )  ,  (231) 

where  j  =  1  or  2  indicates  whether  0og  or  q\og  is  more  accurate,  while  n  is  the  number  of  final 
evaluations  utilized  for  MATLAB  routine  quad/ . 

A  general  review  of  steepest  descent  procedures,  including  an  accurate  prescription  for 
numerically  evaluating  the  paths  of  steepest  descent  out  of  a  saddlepoint,  has  been  presented. 

By  factoring  out  the  value  of  the  integrand  at  the  saddlepoint,  underflow  and  overflow  can  be 
controlled  by  computing  the  logarithm  of  the  specified  quantity.  If  desired,  the  exponential  of 
this  logarithmic  quantity  can  then  be  taken  as  the  last  step  of  computation;  of  course,  this 
exponentiation  could  itself  underflow  or  overflow. 
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APPENDIX 

PROGRAM  LISTING  FOR  QM(a,  b ) 

function  [Q  q  Qlog  qlog  j  n]  =  Qq(M,a,b)  %  Only  for  product  a  b  real. 
%  QM(a,b)  =  int(b,inf)  dx  x  (x/a)/'(M-l)  exp  (-xA2/2-aA2/2)  I  [M-l  ]  (ax)  . 
%  qM(a,b)  =  1  -  QM(a,b)  is  the  complementary  QM  function. 

%  j=l:  Qlog  is  most  accurate;  j  =2 :  qlog  is  most  accurate. 

%  n  is  the  number  of  function  evaluations  in  the  calls  to  quadl . 
global  m  p  r  cs  R  zp  yy 

tolrel=le-6;  %  Relative  Tolerance  on  quadl  integrals. 

R=.2;  %  Threshold  for  SP  &  pole  separation. 

j=i; 

n=0;  %  Number  of  function  evaluations  in  quadl. 

if(b==0)  Q=l;  q=0;  Qlog=0;  qlog=-Inf;  return,  end 
if  (a==0) 

c=b*b/2;  t=l;  s=l; 
for  k=l:M-l 

t=t*c/k;  s=s+t ; 

end 

Q=exp(-c)*s;  q=l-Q;  Qlog=imag_limit (-c+log (s) ) ;  qlog=log(q); 
return 

end 

m=M; 

p=a*b; 

if (abs (imag (p) ) >abs (real (p) ) *le-14)  error ('a*b  is  not  real.*)/  end 
p=real (p) ;  %  Eliminate  imag  round-off  error. 

r=a/b; 

zp=-i*log (r ) ;  %  zp  =  pole  location. 

e=roots ( [p  2*m-2-p*r  -p-2*m*r  p*r] ) ; 

zs=-i*log (e) ;  %  3  SP  (saddlepoint)  locations. 

[d, k] =min (abs (e-r) ) ;  %  [d,k]=min(abs(zs-zp)); 

zs(k)=[];  %  Eliminate  SP  closest  to  pole. 

ang= . 5* (pi-angle (psi2 ( zs ))) ;  %  Retain  the  SP 

al=min (ang (1) ,pi-ang (1) ) ;  %  with  the  lower 

a2=min (ang (2 ) , pi-ang (2 ) ) ;  %  magnitude-slope 

k=l;  if(a2<al)  k=2;  end  %  of  the  POSD 

zs=zs(k);  %  through  the  SP. 

cs=p*cos (zs) +i*m*zs-log (exp (i*zs) -r) ; 
faclog=- . 5* (a*a+b*b) - (m-l) *log (r) -log (2*pi) ; 
xp=real(zp);  yp=imag(zp); 
xs=real (zs) ;  ys=imag(zs);  yy=ys-yp; 
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if (abs (yy) >=R) 

[int, n] =quadl (@f , zs-pi, zs+pi, tolrel) ; 

else 

sq=realsqrt (R*R-yy*yy) ; 

xl=xp-sq;  x2=xp+sq; 

zl=x2-2*pi+i*ys ;  z2=xl+i*ys; 

fm=max (abs ( f ( [ zl  z2]))); 

if  (abs (xs-xp) >=sq)  fm=max ( fm, 1 ) ;  end 

tolabs=fm*tolrel; 

[inti,  nl] =quadl (@f , zl, z2, tolabs) ; 
tp=acos (abs (yy) /R) ; 

t=linspace (-tp, tp, 5 ) ;  %  5  samples  in  angle. 

big=max (abs (g  (t) ) ) ; 

gtolabs=tp*big*tolrel ;  %  Absolute  tolerance  on  quadl  g  integral. 

[ int 2, n2 ] =quadl (@g, -tp, tp, gtolabs) ; 
int=intl+int2 ;  n=nl+n2; 

end 

if (yy<0) 

Qlog=imag_limit (faclog+cs+log (int) ) ; 

Q=exp (Qlog) ;  q=l-Q;  qlog=log(q);  j=l; 

else 

qlog=imag__limit  (faclog+cs  +  log  (int)  +i*pi)  ; 
q=exp (qlog) ;  Q=l-q;  Qlog=log(Q);  j=2; 

end 

if (imag (a) ==0) 

Q=real (Q) ;  q=real (q) ;  Qlog=real (Qlog) ;  qlog=real (qlog) ; 

end 

%  keyboard 


function  w  =  f(z)  %  f(zs)  =  1 

global  m  p  r  cs 

w=exp (p*cos (z) +i*m*z-log  (exp  (i*z) -r) -cs) ; 

function  w  =  g(t) 
global  m  p  r  cs  R  zp  yy 
c=cos(t);  s=sin(t); 
u=R* (s  +  i*c)  ; 

if (yy<0)  u=conj (u) ;  s=-s;  end 
z=zp+u; 

w=R* (c-i*s) . *exp (p*cos (z) +i*m*z-log (exp (i*z) -r) -cs) ; 

function  w  =  psi2(z) 
global  p  r 
e=exp (i*z)  ; 

w=-p*cos (z) -r*e . / (e-r) .A2; 
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